!>  \file rascnv.F90
!!  This file contains the entire Relaxed Arakawa-Schubert convection
!!  parameteriztion

!>\defgroup rascnv_schm Relaxed Arakawa-Schubert Convection Scheme
      module rascnv

      USE machine , ONLY : kind_phys
      implicit none
      public :: rascnv_init, rascnv_run
      private
      logical :: is_initialized = .False.
!
      integer,               parameter :: kp = kind_phys
      integer,               parameter :: nrcmax=32 ! Maximum # of random clouds per 1200s

      integer,               parameter :: idnmax=999
      real (kind=kind_phys), parameter :: delt_c=1800.0_kp/3600.0_kp          &
!     Adjustment time scales in hrs for deep and shallow clouds
!    &,                                   adjts_d=3.0, adjts_s=0.5
!    &,                                   adjts_d=2.5, adjts_s=0.5
     &,                                   adjts_d=2.0_kp, adjts_s=0.5_kp
!
      logical,               parameter :: fix_ncld_hr=.true.

!
      real (kind=kind_phys), parameter :: ZERO=0.0_kp,      HALF=0.5_kp       &
     &,                                   pt25=0.25_kp,     ONE=1.0_kp        &
     &,                                   TWO=2.0_kp,       FOUR=4.0_kp       &
     &,                                   twoo3=two/3.0_kp                    &
     &,                                   FOUR_P2=4.0e2_kp, ONE_M10=1.0e-10_kp&
     &,                                   ONE_M6=1.0e-6_kp, ONE_M5=1.0e-5_kp  &
     &,                                   ONE_M2=1.0e-2_kp, ONE_M1=1.0e-1_kp  &
     &,                                   oneolog10=one/log(10.0_kp)          &
     &,                                   rain_min=1.0e-13_kp                 &
     &,                                   facmb=0.01_kp                       & ! conversion factor from Pa to hPa (or mb)
     &,                                   cmb2pa=100.0_kp                       ! Conversion from hPa to Pa
!
!     real (kind=kind_phys), parameter :: frac=0.5_kp,    crtmsf=0.0_kp     &
      real (kind=kind_phys), parameter :: frac=0.1_kp,    crtmsf=0.0_kp     &
     &,                                   tfrac_max=0.15_kp                 &
     &,                                   rhfacs=0.75_kp, rhfacl=0.75_kp    &
     &,                                   face=5.0_kp,    delx=10000.0_kp   &
     &,                                   ddfac=face*delx*0.001_kp          &
     &,                                   max_neg_bouy=0.15_kp              &
!    &,                                   max_neg_bouy=pt25_kp              &
     &,                                   testmb=0.1_kp, testmbi=one/testmb &
     &,                                   dpd=0.5_kp, rknob=1.0_kp, eknob=1.0_kp

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     logical, parameter :: aw_scal=.false., cumfrc=.true.              &
      logical, parameter :: aw_scal=.true., cumfrc=.true.               &
     &,                     updret=.false., vsmooth=.false.             &
     &,                     wrkfun=.false., crtfun=.true.               &
     &,                     calkbl=.true.,  botop=.true.,  revap=.true. &
     &,                     advcld=.true.,  advups=.false.,advtvd=.true.
!    &,                     advcld=.true.,  advups=.true., advtvd=.false.
!    &,                     advcld=.true.,  advups=.false.,advtvd=.false.


      real(kind=kind_phys), parameter :: TF=233.16_kp, TCR=273.16_kp    &
     &,                                  TCRF=one/(TCR-TF), TCL=2.0_kp

!
!    For pressure gradient force in momentum mixing
!     real (kind=kind_phys), parameter :: pgftop=0.80, pgfbot=0.30       &
!    No pressure gradient force in momentum mixing
      real (kind=kind_phys), parameter :: pgftop=0.0_kp, pgfbot=0.0_kp   &
!     real (kind=kind_phys), parameter :: pgftop=0.55, pgfbot=0.55       &
     &,                                  pgfgrad=(pgfbot-pgftop)*0.001_kp&
     &,                                  cfmax=0.1_kp
!
!     For Tilting Angle Specification
!
      real(kind=kind_phys)       :: REFP(6), REFR(6), TLAC(8), PLAC(8), &
                                    TLBPL(7), drdp(5)
!
      DATA PLAC/100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0/
      DATA TLAC/ 35.0,  25.0,  20.0,  17.5,  15.0,  12.5,  10.0,  7.5/
      DATA REFP/500.0, 300.0, 250.0, 200.0, 150.0, 100.0/
      DATA REFR/  1.0,   2.0,   3.0,   4.0,   6.0,   8.0/
!
      real(kind=kind_phys)       :: AC(16), AD(16)
!
      integer, parameter :: nqrp=500001
      real(kind=kind_phys)       ::  C1XQRP, C2XQRP, TBQRP(NQRP),       &
                                     TBQRA(NQRP), TBQRB(NQRP)
!
      integer, parameter :: nvtp=10001
      real(kind=kind_phys)       ::  C1XVTP, C2XVTP, TBVTP(NVTP)
!
      real(kind=kind_phys)       ::  afc, facdt,                        &
                            grav, cp, alhl, alhf, rgas, rkap, nu, pi,   &
                            t0c,  rv, cvap, cliq, csol, ttp, eps, epsm1,&
!
                            ONEBG,   GRAVCON,  onebcp, GRAVFAC, ELOCP,  &
                            ELFOCP,  oneoalhl, CMPOR,  picon,   zfac,   &
                            deg2rad, PIINV,    testmboalhl,             &
                            rvi,     facw,     faci, hsub, tmix, DEN

      contains
 
! -----------------------------------------------------------------------
! CCPP entry points for gfdl cloud microphysics
! -----------------------------------------------------------------------

!>\ingroup rascnv_schm
!>\brief The subroutine initializes rascnv
!!
!> \section arg_table_rascnv_init Argument Table
!! \htmlinclude rascnv_init.html
!!
      subroutine rascnv_init(me, dt,   con_g,    con_cp,   con_rd,      &
                             con_rv,   con_hvap, con_hfus, con_fvirt,   &
                             con_t0c,  con_ttp,  con_cvap, con_cliq,    &
                             con_csol, con_eps,  con_epsm1,             &
                             errmsg,   errflg)
!
      Implicit none
!
      integer,              intent(in)  :: me
      real(kind=kind_phys), intent(in)  :: dt,                          &
                      con_g,    con_cp,    con_rd,  con_rv,   con_hvap, &
                      con_hfus, con_fvirt, con_t0c, con_cvap, con_cliq, &
                      con_csol, con_ttp,   con_eps, con_epsm1

      character(len=*),     intent(out) :: errmsg
      integer,              intent(out) :: errflg
!
      real(kind=kind_phys), parameter ::  actp=1.7_kp, facm=1.00_kp
!
      real(kind=kind_phys) :: PH(15), A(15)
!
      DATA PH/150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0    &
     &,       550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0/
!
       DATA A/ 1.6851, 1.1686, 0.7663, 0.5255, 0.4100, 0.3677           &
     &,        0.3151, 0.2216, 0.1521, 0.1082, 0.0750, 0.0664           &
     &,        0.0553, 0.0445, 0.0633/
!
      real(kind=kind_phys) tem,  actop, tem1, tem2
      integer              i, l
!
! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0
      if (is_initialized) return
!                           set critical workfunction arrays
      ACTOP = ACTP*FACM
      DO L=1,15
        A(L) = A(L)*FACM
      ENDDO
      DO L=2,15
        TEM   = one / (PH(L) - PH(L-1))
        AC(L) = (PH(L)*A(L-1) - PH(L-1)*A(L)) * TEM
        AD(L) = (A(L) - A(L-1)) * TEM
      ENDDO
      AC(1)  = ACTOP
      AC(16) = A(15)
      AD(1)  = zero
      AD(16) = zero
!
      CALL SETQRP
      CALL SETVTP
!
      do i=1,7
        tlbpl(i) = (tlac(i)-tlac(i+1)) / (plac(i)-plac(i+1))
      enddo
      do i=1,5
        drdp(i)  = (REFR(i+1)-REFR(i)) / (REFP(i+1)-REFP(i))
      enddo
!
!     VTP = 36.34*SQRT(1.2)* (0.001)**0.1364
!
      AFC = -(1.01097e-4_kp*DT)*(3600.0_kp/DT)**0.57777778_kp
!
      if (fix_ncld_hr) then
        facdt = delt_c / dt
      else
        facdt = one / 3600.0_kp
      endif
!
      grav = con_g     ; cp    = con_cp   ; alhl = con_hvap
      alhf = con_hfus  ; rgas  = con_rd
      nu   = con_FVirt ; t0c   = con_t0c
      rv   = con_rv    ; cvap  = con_cvap
      cliq = con_cliq  ; csol  = con_csol ; ttp  = con_ttp
      eps  = con_eps   ; epsm1 = con_epsm1
!
      pi       = four*atan(one) ; PIINV   = one/PI
      ONEBG    = ONE / GRAV     ; GRAVCON = cmb2pa * ONEBG
      onebcp   = one / cp       ; GRAVFAC = GRAV / CMB2PA
      rkap     = rgas * onebcp  ; deg2rad = pi/180.0_kp
      ELOCP    = ALHL * onebcp  ; ELFOCP  = (ALHL+ALHF) * onebcp
      oneoalhl = one/alhl       ; CMPOR   = CMB2PA / RGAS
      picon    = half*pi*onebg  ; zfac    = 0.28888889e-4_kp * ONEBG
      testmboalhl = testmb/alhl
!
      rvi  = one / rv       ; facw = CVAP - CLIQ 
      faci = CVAP - CSOL    ; hsub = alhl + alhf
      tmix = TTP  - 20.0_kp ; DEN  = one / (TTP-TMIX)
!

      if (me == 0) write(0,*) ' NO DOWNDRAFT FOR CLOUD TYPES'           &
     &,                ' DETRAINING AT NORMALIZED PRESSURE ABOVE ',DPD
!
      is_initialized = .true.

!
      end subroutine rascnv_init
!
!!
!!===================================================================== !
!! rascnv_run:                                                          !
!!                                                                      !
!! program history log:                                                 !
!!   Oct  2019  -- shrinivas moorthi                                    !
!!                                                                      !
!!                                                                      !
!! ====================  defination of variables  ====================
!! !
!!                                                                      !
!! inputs:                                                       size
!! !
!!    im       - integer, horiz dimension and num of used pts      1    !
!!    k        - integer, vertical dimension                       1    !
!!    dt       - real, time step in seconds                        1    !
!!    dtf      - real, dynamics time step in seconds               1    !
!!    rannum   - real, array holding random numbers between 0 an 1 (im,nrcm)  !
!!    tin      - real, input temperature (K)
!!    qin      - real, input specific humidity (kg/kg)
!!    uin      - real, input zonal wind component
!!    vin      - real, input meridional wind component
!!    ccin     - real, input condensates+tracers 
!!    fscav    - real 
!!    prsi     - real, layer interface pressure
!!    prsl     - real, layer mid pressure
!!    prsik    - real, layer interface Exner function
!!    prslk    - real, layer mid Exner function
!!    phil     - real, layer mid geopotential height
!!    phii     - real, layer interface geopotential height
!!    kpbl     - integer pbl top index
!!    cdrag    - real, drag coefficient
!!    rainc    - real, convectinve rain (m/sec)
!!    kbot     - integer, cloud bottom index
!!    ktop     - integer, cloud top index
!!    knv      - integer, 0 - no convvection; 1 - convection
!!    ddvel    - downdraft induced surface wind
!!    flipv    - logical, true if input data from bottom to top
!!    me       - integer, current pe number
!!    area     - real, grid area
!!    ccwf     - real, multiplication factor for critical workfunction
!!    nrcm     - integer, number of random numbers at each grid point
!!    rhc      - real, critical relative humidity
!!    ud_mf    - real, updraft mass flux
!!    dd_mf    - real, downdraft mass flux
!!    dt_mf    - real, detrained mass flux
!!    qw0      - real, min cloud water before autoconversion
!!    qi0      - real, min cloud ice before autoconversion
!!    dlqfac   - real,fraction of condensated detrained in layers
!!    kdt      - integer, current teime step
!!    revap    - logial,  when true reevaporate falling rain/snow
!!    qlcn     - real
!!    qicn     - real
!!    w_upi    - real
!!    cf_upi   - real
!!    cnv_mfd  - real
!!    cnv_dqldt- real
!!    clcn     - real
!!    cnv_fice - real
!!    cnv_ndrop- real
!!    cnv_nice - real
!!    mp_phys  - integer, microphysics option
!!    mp_phys_mg - integer, flag for MG microphysics option
!!    trcmin   - real, floor value for tracers
!!    ntk      - integer, index representing TKE in the tracer array
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!>\ingroup rascnv_schm
!! \section arg_table_rascnv_run Argument Table
!! \htmlinclude rascnv_run.html
!!
      subroutine rascnv_run(IM,     k,      itc, ntc, ntr,   dt,  dtf   &
     &,                     ccwf,   area,   dxmin, dxinv                &
     &,                     psauras, prauras, wminras, dlqf, flipv      &
     &,                     me,     rannum, nrcm,  mp_phys, mp_phys_mg  &
     &,                     ntk,    kdt,                            rhc &
     &,                     tin,    qin,    uin,   vin,     ccin, fscav &
     &,                     prsi,   prsl,   prsik, prslk,   phil, phii  &
     &,                     KPBL,   CDRAG,  RAINC, kbot,    ktop, kcnv  &
     &,                     DDVEL,  ud_mf,  dd_mf, dt_mf                &
     &,                     QLCN,   QICN,   w_upi,  cf_upi, CNV_MFD     &
     &,                     CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE  &
     &,                     errmsg,  errflg)
!
!*********************************************************************
!*********************************************************************
!************         Relaxed Arakawa-Schubert      ******************
!************             Parameterization          ******************
!************          Plug Compatible Driver       ******************
!************               23 May 2002             ******************
!************                                       ******************
!************               Developed By            ******************
!************                                       ******************
!************             Shrinivas Moorthi         ******************
!************                                       ******************
!************                  EMC/NCEP             ******************
!*********************************************************************
!*********************************************************************
!
!
      Implicit none
!
!      input
!
      logical, intent(in) :: flipv
!
      integer, intent(in) :: im, k, itc, ntc, ntr, me, nrcm, ntk, kdt   &
     &,                      mp_phys, mp_phys_mg
      integer, dimension(:), intent(out)   :: kbot, ktop
      integer, dimension(:), intent(inout) :: kcnv
      integer, dimension(:), intent(in)    :: kpbl
!
      real(kind=kind_phys), intent(in)        :: dxmin, dxinv, ccwf(:)  &
     &,                                          psauras(:), prauras(:) &
     &,                                          wminras(:), dlqf(:)
!
      real(kind=kind_phys), dimension(:,:), intent(in) :: prsi, prsik, phii

      real(kind=kind_phys), dimension(:,:), intent(inout) :: tin, qin,  uin, vin
      real(kind=kind_phys), dimension(:,:), intent(in) :: prsl, prslk, phil      &
     &,                                                   rhc
      real(kind=kind_phys), dimension(:,:), intent(out) :: ud_mf, dd_mf, dt_mf
      real(kind=kind_phys), dimension(:,:), intent(inout) :: qlcn, qicn, w_upi   &
     &,                                                   cnv_mfd                &
     &,                                                   cnv_dqldt, clcn        &
     &,                                                   cnv_fice, cnv_ndrop    &
     &,                                                   cnv_nice, cf_upi
      real(kind=kind_phys), dimension(:)  , intent(in)  :: area,  cdrag
      real(kind=kind_phys), dimension(:)  , intent(out) :: rainc, ddvel
      real(kind=kind_phys), dimension(:,:), intent(in)  :: rannum
      real(kind=kind_phys), intent(inout) :: ccin(:,:,:)
      real(kind=kind_phys), intent(in)    :: dt, dtf
!
!     Added for aerosol scavenging for GOCART
!
      real(kind=kind_phys), intent(in)  :: fscav(:)

!    &,                                   ctei_r(im), ctei_rm
      character(len=*),     intent(out) :: errmsg
      integer,              intent(out) :: errflg
!
!     locals
!
      real(kind=kind_phys)                 :: trcmin(ntr+2)
      real(kind=kind_phys), dimension(k)   :: toi,    qoi, tcu, qcu     &
     &,                                       pcu,    clw, cli, qii, qli&
     &,                                       phi_l,  prsm,psjm         &
     &,                                               alfind, rhc_l     &
!    &,                                       alfinq, alfind, rhc_l     &
     &,                                       qoi_l, qli_l, qii_l
      real(kind=kind_phys), dimension(k+1) :: prs, psj, phi_h, flx, flxd


      integer, dimension(100)         :: ic
      real(kind=kind_phys), parameter :: clwmin=1.0e-10_kp
!
      real(kind=kind_phys), allocatable ::  ALFINT(:,:), uvi(:,:)       &
     &,                                     trcfac(:,:), rcu(:,:)
      real(kind=kind_phys)                  dtvd(2,4)
!    &,                                     DPI(K)
      real(kind=kind_phys) CFAC, TEM,  sgc, ccwfac, tem1, tem2, rain    &
     &,                    wfnc,tla,pl,qiid,qlid, c0, c0i, dlq_fac, sumq&
     &,                    rainp
!     integer                          :: nrcmax    ! Maximum # of random clouds per 1200s
!
      Integer              KCR,  KFX, NCMX, NC,  KTEM, I,   ii, Lm1, l  &
     &,                    ntrc,      ll,   km1, kp1,  ipt, lv, KBL, n  &
     &,                    KRMIN, KRMAX, KFMAX, kblmx, irnd,ib          &
     &,                    kblmn, ksfc, ncrnd
      real(kind=kind_phys) sgcs(k)
!
!  Scavenging related parameters
!
      real                fscav_(ntr+2)  ! Fraction scavenged per km
!
      fscav_ = -999.0_kp                 ! By default no scavenging
      if (itc > 0 .and. ntc > 0) then
        n = itc + ntc - 1
        if (n <= ntr + 2) then
          fscav_(itc:n) = fscav
        else
          errmsg = 'Error in rascnv_run: test ntr >= itc + ntc - 3 FAILED'
          errflg = 1
          return
        end if
      end if
      trcmin = -99999.0_kp
      if (ntk-2 > 0) trcmin(ntk-2) = 1.0e-4_kp

!> - Initialize CCPP error handling variables

      errmsg = ''
      errflg = 0
!
      km1 = k - 1
      kp1 = k + 1
      if (flipv) then
        ksfc = 1
      else
        ksfc = kp1
      endif
!
      ntrc = ntr
      IF (CUMFRC) THEN
        ntrc = ntrc + 2
      ENDIF
      if (ntrc > 0) then
        if (.not. allocated(trcfac)) allocate (trcfac(k,ntrc))
        if (.not. allocated(uvi))    allocate (uvi(k,ntrc))
        if (.not. allocated(rcu))    allocate (rcu(k,ntrc))
        do n=1, ntrc
          do l=1,k
            trcfac(l,n) = one         !  For other tracers
            rcu(l,n)    = zero
          enddo
        enddo
      endif
!
!!!!! initialization for microphysics ACheng
      if(mp_phys == mp_phys_mg) then
        do l=1,K
          do i=1,im
            QLCN(i,l)      = zero
            QICN(i,l)      = zero
            w_upi(i,l)     = zero
            cf_upi(i,l)    = zero
            CNV_MFD(i,l)   = zero
!           CNV_PRC3(i,l)  = zero
            CNV_DQLDT(i,l) = zero
            CLCN(i,l)      = zero
            CNV_FICE(i,l)  = zero
            CNV_NDROP(i,l) = zero
            CNV_NICE(i,l)  = zero
          enddo
        enddo
      endif
!
      if (.not. allocated(alfint)) allocate(alfint(k,ntrc+4))
!
!     call set_ras_afc(dt)
!     AFC = -(1.04E-4*DT)*(3600./DT)**0.578
!     AFC = -(1.01097E-4*DT)*(3600./DT)**0.57777778
!
      do l=1,k
        do i=1,im
          ud_mf(i,l) = zero
          dd_mf(i,l) = zero
          dt_mf(i,l) = zero
        enddo
      enddo
      DO IPT=1,IM

        tem1    = max(zero, min(one, (log(area(ipt)) - dxmin) * dxinv))
        tem2    = one - tem1
        ccwfac  = ccwf(1)*tem1 + ccwf(2)*tem2
        dlq_fac = dlqf(1)*tem1 + dlqf(2)*tem2
        tem     = one  + dlq_fac
        c0i     = (psauras(1)*tem1 + psauras(2)*tem2) * tem
        c0      = (prauras(1)*tem1 + prauras(2)*tem2) * tem
        if (ccwfac == zero) ccwfac = half
!
!       ctei = .false.
!       if (ctei_r(ipt) > ctei_rm) ctei = .true.
!
!     Compute NCRND  : 
!                      if flipv is true, then input variables are from bottom
!                      to top while RAS goes top to bottom
!
        tem = one / prsi(ipt,ksfc)

        KRMIN = 1
        KRMAX = km1
        KFMAX = KRMAX
        kblmx = 1
        kblmn = 1
        sgcs(k) = one
        DO L=1,KM1
          ll = l
          if (flipv) ll = kp1 -l ! Input variables are bottom to top!
          SGC = prsl(ipt,ll) * tem
          sgcs(l) = sgc
          IF (SGC <= 0.050_kp) KRMIN = L
!         IF (SGC <= 0.700_kp) KRMAX = L
!         IF (SGC <= 0.800_kp) KRMAX = L
          IF (SGC <= 0.760_kp) KRMAX = L
!         IF (SGC <= 0.930_kp) KFMAX = L
          IF (SGC <= 0.970_kp) KFMAX = L    ! Commented on 20060202
!         IF (SGC <= 0.700_kp) kblmx = L    ! Commented on 20101015
          IF (SGC <= 0.600_kp) kblmx = L    ! 
!         IF (SGC <= 0.650_kp) kblmx = L    ! Commented on 20060202
          IF (SGC <= 0.980_kp) kblmn = L    ! 
        ENDDO
        krmin = max(krmin,2)

!
        if (fix_ncld_hr) then
!!!       NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/1200) + 0.50001
          NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/1800) + 0.10001_kp
!         NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/1200) + 0.10001
!         NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/900) + 0.50001
!         NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/600) + 0.50001
!         NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * (DTF/360) + 0.50001
!    &                                         + 0.50001
!         NCRND = min(nrcmax, (KRMAX-KRMIN+1)) * min(1.0,DTF/360) + 0.1
        else
          NCRND = min(nrcmax, (KRMAX-KRMIN+1))
        endif
        NCRND   = min(nrcm,max(NCRND, 1))
!
        KCR     = MIN(K,KRMAX)
        KTEM    = MIN(K,KFMAX)
        KFX     = KTEM - KCR

        IF (KFX > 0) THEN
          IF (BOTOP) THEN
            DO NC=1,KFX
              IC(NC) = KTEM + 1 - NC
            ENDDO
          ELSE
            DO NC=KFX,1,-1
              IC(NC) = KTEM + 1 - NC
            ENDDO
          ENDIF
        ENDIF
!
        NCMX  = KFX + NCRND
        IF (NCRND > 0) THEN
          DO I=1,NCRND
            II = mod(i-1,nrcm) + 1
            IRND = (RANNUM(ipt,II)-0.0005_kp)*(KCR-KRMIN+1)
            IC(KFX+I) = IRND + KRMIN
          ENDDO
        ENDIF
!
        do l=1,k
          CLW(l)     = zero
          CLI(l)     = zero
                                 ! to be zero i.e. no environmental condensate!!!
          QII(l)     = zero
          QLI(l)     = zero
!                          Initialize heating, drying, cloudiness etc.
          tcu(l)     = zero
          qcu(l)     = zero
          pcu(l)     = zero
          flx(l)     = zero
          flxd(l)    = zero
          do n=1,ntrc
            rcu(l,n) = zero
          enddo
        enddo
        flx(kp1)  = zero
        flxd(kp1) = zero
        rain      = zero
!
        if (flipv) then   ! Input variables are bottom to top!
          do l=1,k
            ll = kp1 - l
                          ! Transfer input prognostic data into local variable
            toi(l)   = tin(ipt,ll)
            qoi(l)   = qin(ipt,ll)

            PRSM(L)  = prsl(ipt,ll) * facmb
            PSJM(L)  = prslk(ipt,ll)
            phi_l(L) = phil(ipt,ll)
            rhc_l(L) = rhc(ipt,ll)
!
            if (ntrc > ntr) then               ! CUMFRC is true
              uvi(l,ntr+1) = uin(ipt,ll)
              uvi(l,ntr+2) = vin(ipt,ll)
            endif
!
            if (ntr > 0) then                  ! tracers such as O3, dust etc
              do n=1,ntr
                uvi(l,n) = ccin(ipt,ll,n+2)
                if (abs(uvi(l,n)) < 1.0e-20_kp) uvi(l,n) = zero
              enddo
            endif
          enddo
          do l=1,kp1
            ll = kp1 + 1 - l      ! Input variables are bottom to top!
            PRS(LL)   = prsi(ipt,L) * facmb
            PSJ(LL)   = prsik(ipt,L)
            phi_h(LL) = phii(ipt,L)
          enddo
!
          if (ccin(ipt,1,2) <= -998.0_kp) then  ! input ice/water are together 
            do l=1,k
              ll = kp1 -l
              tem = ccin(ipt,ll,1)                                      &
     &            * MAX(ZERO, MIN(ONE, (TCR-toi(L))*TCRF))
              ccin(ipt,ll,2) = ccin(ipt,ll,1) - tem
              ccin(ipt,ll,1) = tem
            enddo
          endif
          if (advcld) then
            do l=1,k
              ll = kp1 -l ! Input variables are bottom to top!
              QII(L) = ccin(ipt,ll,1)
              QLI(L) = ccin(ipt,ll,2)
            enddo
          endif
          KBL  = MAX(MIN(k, kp1-KPBL(ipt)), k/2)
!
        else              ! Input variables are top to bottom!

          do l=1,k
                          ! Transfer input prognostic data into local variable
            toi(l)   = tin(ipt,l)
            qoi(l)   = qin(ipt,l)

            PRSM(L)  = prsl(ipt, L) * facmb
            PSJM(L)  = prslk(ipt,L)
            phi_l(L) = phil(ipt,L)
            rhc_l(L) = rhc(ipt,L)
!
            if (ntrc > ntr) then               ! CUMFRC is true
              uvi(l,ntr+1) = uin(ipt,l)
              uvi(l,ntr+2) = vin(ipt,l)
            endif
!
            if (ntr > 0) then                  ! tracers such as O3, dust etc
              do n=1,ntr
                uvi(l,n) = ccin(ipt,l,n+2)
                if (abs(uvi(l,n)) < 1.0e-20_kp) uvi(l,n) = zero
              enddo
            endif
          enddo
          DO L=1,kp1
            PRS(L)   = prsi(ipt,L) * facmb
            PSJ(L)   = prsik(ipt,L)
            phi_h(L) = phii(ipt,L)
          ENDDO
!
          if (ccin(ipt,1,2) <= -998.0_kp) then  ! input ice/water are together
            do l=1,k
              tem = ccin(ipt,l,1)                                       &
     &            * MAX(ZERO, MIN(ONE, (TCR-toi(L))*TCRF))
              ccin(ipt,l,2) = ccin(ipt,l,1) - tem
              ccin(ipt,l,1) = tem
            enddo
          endif
          if (advcld) then
            do l=1,k
              QII(L) = ccin(ipt,l,1)
              QLI(L) = ccin(ipt,l,2)
            enddo
          endif
!
          KBL  = KPBL(ipt)
!
        endif      ! end of if (flipv) then
!
!       do l=k,kctop(1),-1
!!        DPI(L)  = 1.0 / (PRS(L+1) - PRS(L))
!       enddo
!
!     print *,' ipt=',ipt

        if (advups) then               ! For first order upstream for updraft
          alfint(:,:) = one
        elseif (advtvd) then           ! TVD flux limiter scheme for updraft
!         alfint(:,:) = one
          alfint(:,:) = half
          l   = krmin
          lm1 = l - 1
          dtvd(1,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)        &
     &              + alhl*(qoi(l)-qoi(lm1))
          dtvd(1,2) = qoi(l) - qoi(lm1)
          dtvd(1,3) = qli(l) - qli(lm1)
          dtvd(1,4) = qii(l) - qii(lm1)
          do l=krmin+1,k
            lm1 = l - 1

!     write(0,*)' toi=',toi(l),toi(lm1),' phi_l=',phi_l(l),phi_l(lm1)
!    &,' qoi=',qoi(l),qoi(lm1),' cp=',cp,' alhl=',alhl

            dtvd(2,1)   = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)    &
     &                  + alhl*(qoi(l)-qoi(lm1))

!     write(0,*)' l=',l,' dtvd=',dtvd(:,1)

            if (abs(dtvd(2,1)) > 1.0e-10_kp) then
              tem1        = dtvd(1,1) / dtvd(2,1)
              tem2        = abs(tem1)
              alfint(l,1) = one - half*(tem1 + tem2)/(one + tem2)   ! for h
            endif

!     write(0,*)' alfint=',alfint(l,1),' l=',l,' ipt=',ipt

            dtvd(1,1) = dtvd(2,1)
!
            dtvd(2,2) = qoi(l) - qoi(lm1)

!     write(0,*)' l=',l,' dtvd2=',dtvd(:,2)

            if (abs(dtvd(2,2)) > 1.0e-10_kp) then
              tem1        = dtvd(1,2) / dtvd(2,2)
              tem2        = abs(tem1)
              alfint(l,2) = one - half*(tem1 + tem2)/(one + tem2)   ! for q
            endif
            dtvd(1,2) = dtvd(2,2)
!
            dtvd(2,3) = qli(l) - qli(lm1)

!     write(0,*)' l=',l,' dtvd3=',dtvd(:,3)

            if (abs(dtvd(2,3)) > 1.0e-10_kp) then
              tem1        = dtvd(1,3) / dtvd(2,3)
              tem2        = abs(tem1)
              alfint(l,3) = one - half*(tem1 + tem2)/(one + tem2)   ! for ql
            endif
            dtvd(1,3) = dtvd(2,3)
!
            dtvd(2,4) = qii(l) - qii(lm1)

!     write(0,*)' l=',l,' dtvd4=',dtvd(:,4)

            if (abs(dtvd(2,4)) > 1.0e-10_kp) then
              tem1        = dtvd(1,4) / dtvd(2,4)
              tem2        = abs(tem1)
              alfint(l,4) = one - half*(tem1 + tem2)/(one + tem2)   ! for qi
            endif
            dtvd(1,4) = dtvd(2,4)
          enddo
!
          if (ntrc > 0) then
            do n=1,ntrc
              l = krmin
              dtvd(1,1) = uvi(l,n) - uvi(l-1,n)
              do l=krmin+1,k
                dtvd(2,1) = uvi(l,n) - uvi(l-1,n)

!     write(0,*)' l=',l,' dtvdn=',dtvd(:,1),' n=',n,' l=',l

                if (abs(dtvd(2,1)) > 1.0e-10_kp) then
                  tem1          = dtvd(1,1) / dtvd(2,1)
                  tem2          = abs(tem1)
                  alfint(l,n+4) = one - half*(tem1 + tem2)/(one + tem2) ! for tracers
                endif
                dtvd(1,1) = dtvd(2,1)
              enddo
            enddo
          endif
        else
          alfint(:,:) = half              ! For second order scheme
        endif
        alfind(:)   = half
!
!     write(0,*)' after alfint for ipt=',ipt

! Resolution dependent press grad correction momentum mixing

        if (CUMFRC) then
          do l=krmin,k
            tem = one - max(pgfbot, min(pgftop, pgftop+pgfgrad*prsm(l)))
            trcfac(l,ntr+1) = tem
            trcfac(l,ntr+2) = tem
          enddo
        endif
!
!       if (calkbl) kbl = k

        if (calkbl) then
          kbl = kblmn
        else
          kbl = min(kbl, kblmn)
        endif
!
        DO NC=1,NCMX     ! multi cloud loop
!
          IB = IC(NC)    ! cloud top level index
          if (ib > kbl-1) cycle
!
!****************************************************************************
!         if (advtvd) then           ! TVD flux limiter scheme for updraft
!           l   = ib
!           lm1 = l - 1
!           dtvd(1,1) = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)
!    &                + alhl*(qoi(l)-qoi(lm1))
!           dtvd(1,2) = qoi(l) - qoi(lm1)
!           dtvd(1,3) = qli(l) - qli(lm1)
!           dtvd(1,4) = qii(l) - qii(lm1)
!           do l=ib+1,k
!             lm1 = l - 1
!             dtvd(2,1)   = cp*(toi(l)-toi(lm1)) + phi_l(l)-phi_l(lm1)
!    &                    + alhl*(qoi(l)-qoi(lm1))
!             if (abs(dtvd(2,1)) > 1.0e-10) then
!               tem1        = dtvd(1,1) / dtvd(2,1)
!               tem2        = abs(tem1)
!               alfint(l,1) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for h
!             endif
!             dtvd(1,1)   = dtvd(2,1)
!
!             dtvd(2,2)   = qoi(l) - qoi(lm1)
!             if (abs(dtvd(2,2)) > 1.0e-10) then
!               tem1        = dtvd(1,2) / dtvd(2,2)
!               tem2        = abs(tem1)
!               alfint(l,2) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for q
!             endif
!             dtvd(1,2)   = dtvd(2,2)
!
!             dtvd(2,3)   = qli(l) - qli(lm1)
!             if (abs(dtvd(2,3)) > 1.0e-10) then
!               tem1        = dtvd(1,3) / dtvd(2,3)
!               tem2        = abs(tem1)
!               alfint(l,3) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for ql
!             endif
!             dtvd(1,3)   = dtvd(2,3)
!
!             dtvd(2,4)   = qii(l) - qii(lm1)
!             if (abs(dtvd(2,4)) > 1.0e-10) then
!               tem1        = dtvd(1,4) / dtvd(2,4)
!               tem2        = abs(tem1)
!               alfint(l,4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2)   ! for qi
!             endif
!             dtvd(1,4)   = dtvd(2,4)
!           enddo
!
!           if (ntrc > 0) then
!             do n=1,ntrc
!               l = ib
!               dtvd(1,1)   = uvi(l,n) - uvi(l-1,n)
!               do l=ib+1,k
!                 dtvd(2,1)     = uvi(l,n) - uvi(l-1,n)
!                 if (abs(dtvd(2,1)) > 1.0e-10) then
!                   tem1        = dtvd(1,1) / dtvd(2,1)
!                   tem2          = abs(tem1)
!                   alfint(l,n+4) = 1.0 - 0.5*(tem1 + tem2)/(1.0 + tem2) ! for tracers
!                 endif
!                 dtvd(1,1)     = dtvd(2,1)
!               enddo
!             enddo
!           endif
!         endif
!****************************************************************************
!
          WFNC = zero
          do L=IB,KP1
            FLX(L)  = zero
            FLXD(L) = zero
          enddo
!
          TLA  = -10.0_kp
!
          qiid = qii(ib)         ! cloud top level ice before convection
          qlid = qli(ib)         ! cloud top level water before convection
!
          rainp = rain

          CALL CLOUD(K,  KP1, IB, ntrc, kblmx, kblmn                    &
     &,              FRAC,  MAX_NEG_BOUY, vsmooth, aw_scal              &
     &,              REVAP, WRKFUN, CALKBL, CRTFUN                      &
     &,              DT, KDT, TLA, DPD                                  &
     &,              ALFINT, rhfacl, rhfacs, area(ipt)                  &
     &,              ccwfac, CDRAG(ipt), trcfac                         &
     &,              alfind, rhc_l, phi_l, phi_h, PRS, PRSM,sgcs        &
     &,              TOI, QOI, UVI, QLI, QII, KBL, DDVEL(ipt)           &
     &,              TCU, QCU, RCU, PCU, FLX, FLXD, RAIN, WFNC, fscav_  &
     &,              trcmin, ntk-2, c0, wminras(1), c0i, wminras(2)     &
     &,              dlq_fac)
!    &,              ctei)

!
          if (flipv) then
            do L=IB,K
              ll = kp1 -l    ! Input variables are bottom to top!
              ud_mf(ipt,ll) = ud_mf(ipt,ll) + flx(l+1)
              dd_mf(ipt,ll) = dd_mf(ipt,ll) + flxd(l+1)
            enddo
            ll  = kp1 - ib
            dt_mf(ipt,ll) = dt_mf(ipt,ll) + flx(ib)

            if (mp_phys == mp_phys_mg) then !         Anning Cheng for microphysics 11/14/2015

              CNV_MFD(ipt,ll)   = CNV_MFD(ipt,ll)   + flx(ib)/dt

!!            CNV_DQLDT(ipt,ll) = CNV_DQLDT(ipt,ll)
!!   &                          + max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt
              CNV_DQLDT(ipt,ll) = CNV_DQLDT(ipt,ll) + flx(ib)*          &
     &                            max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt
!!   &                                max(0.,(QLI(ib)+QII(ib)))/dt/3.
              if(flx(ib)<0) write(*,*)"AAA666", flx(ib),QLI(ib),QII(ib) &
     &                                       ,ipt,ll
            endif

          else

            do L=IB,K
              ud_mf(ipt,l) = ud_mf(ipt,l) + flx(l+1)
              dd_mf(ipt,l) = dd_mf(ipt,l) + flxd(l+1)
            enddo
            dt_mf(ipt,ib) = dt_mf(ipt,ib) + flx(ib)

            if (mp_phys == mp_phys_mg) then !         Anning Cheng for microphysics 11/14/2015
              CNV_MFD(ipt,ib)   = CNV_MFD(ipt,ib)   + flx(ib)/dt
!!            CNV_DQLDT(ipt,ib) = CNV_DQLDT(ipt,ib)
!!   &                          + max(0.,(QLI(ib)+QII(ib)-qiid-qlid))/dt
              CNV_DQLDT(ipt,ib) = CNV_DQLDT(ipt,ib) + flx(ib)*          &
     &                            max(zero,(QLI(ib)+QII(ib)-qiid-qlid))/dt
!!   &                                max(0.,(QLI(ib)+QII(ib)))/dt/3.
              if(flx(ib)<0) write(*,*)"AAA666", flx(ib),QLI(ib),QII(ib) &
     &                                       ,ipt,ib
            endif

          endif
! 
!
!   Warning!!!!
!   ------------
!   By doing the following, CLOUD does not contain environmental
!   condensate!
!
          if (.not. advcld) then
            do l=1,K
              clw(l) = clw(l) + QLI(L)
              cli(l) = cli(l) + QII(L)
              QLI(L) = zero
              QII(L) = zero
            enddo
          endif
!
        ENDDO                        ! End of the NC loop!
!
        RAINC(ipt) = rain * 0.001_kp    ! Output rain is in meters
        if (rainc(ipt) < rain_min) rainc(ipt) = zero

        ktop(ipt) = kp1
        kbot(ipt) = 0

        kcnv(ipt) = 0


        do l=k,1,-1
!         qli(l) = max(qli(l), zero)
!         qii(l) = max(qii(l), zero)
!         clw(i) = max(clw(i), zero)
!         cli(i) = max(cli(i), zero)

          if (sgcs(l) < 0.93_kp .and. abs(tcu(l)) > one_m10) then
!         if (sgcs(l) < 0.90_kp .and. tcu(l) .ne. zero) then
!         if (sgcs(l) < 0.85_kp .and. tcu(l) .ne. zero) then
             kcnv(ipt) = 1
          endif
!  New test for convective clouds ! added in 08/21/96
          if (clw(l)+cli(l) > zero .OR.                                 &
     &        qli(l)+qii(l) > clwmin) ktop(ipt) = l
        enddo
        do l=1,km1
          if (clw(l)+cli(l) > zero .OR.                                 &
     &        qli(l)+qii(l) > clwmin) kbot(ipt) = l
        enddo
!
        if (flipv) then
          do l=1,k
            ll = kp1 - l
            tin(ipt,ll) = toi(l)                  ! Temperature
            qin(ipt,ll) = qoi(l)                  ! Specific humidity
            uin(ipt,ll) = uvi(l,ntr+1)            ! U momentum
            vin(ipt,ll) = uvi(l,ntr+2)            ! V momentum

!!        for 2M microphysics, always output these variables
            if (mp_phys == mp_phys_mg) then
              if (advcld) then
                QLCN(ipt,ll)     = max(qli(l)-ccin(ipt,ll,2), zero)
                QICN(ipt,ll)     = max(qii(l)-ccin(ipt,ll,1), zero)
                CNV_FICE(ipt,ll) = QICN(ipt,ll)                         &
     &                           / max(1.0e-10_kp,QLCN(ipt,ll)+QICN(ipt,ll))
              else
                QLCN(ipt,ll)     = qli(l)
                QICN(ipt,ll)     = qii(l)
                CNV_FICE(ipt,ll) = qii(l)/max(1.0e-10_kp,qii(l)+qli(l))
              endif
              cf_upi(ipt,ll)   = max(zero,min(0.02_kp*log(one+           &
     &                             500.0_kp*ud_mf(ipt,ll)/dt), cfmax))
!    &                             500*ud_mf(ipt,ll)/dt), 0.60))
              CLCN(ipt,ll)     = cf_upi(ipt,ll)  !downdraft is below updraft
              w_upi(ipt,ll)    = ud_mf(ipt,ll)*toi(l)*rgas /            &
     &                      (dt*max(cf_upi(ipt,ll),1.0e-12_kp)*prsl(ipt,ll))
            endif

            if (ntr > 0) then
              do n=1,ntr
                ccin(ipt,ll,n+2) = uvi(l,n)           ! Tracers
              enddo
            endif
          enddo
          if (advcld) then
            do l=1,k
              ll  = kp1 - l
              ccin(ipt,ll,1) = qii(l)          ! Cloud ice
              ccin(ipt,ll,2) = qli(l)          ! Cloud water
            enddo
          else
            do l=1,k
              ll  = kp1 - l
              ccin(ipt,ll,1) = ccin(ipt,ll,1) + cli(l)
              ccin(ipt,ll,2) = ccin(ipt,ll,2) + clw(l)
            enddo
          endif
!
          ktop(ipt) = kp1 - ktop(ipt)
          kbot(ipt) = kp1 - kbot(ipt)
!
        else

          do l=1,k
            tin(ipt,l) = toi(l)                   ! Temperature
            qin(ipt,l) = qoi(l)                   ! Specific humidity
            uin(ipt,l) = uvi(l,ntr+1)             ! U momentum
            vin(ipt,l) = uvi(l,ntr+2)             ! V momentum

!!        for 2M microphysics, always output these variables
            if (mp_phys == mp_phys_mg) then
              if (advcld) then
                QLCN(ipt,l)     = max(qli(l)-ccin(ipt,l,2), zero)
                QICN(ipt,l)     = max(qii(l)-ccin(ipt,l,1), zero)
                CNV_FICE(ipt,l) = QICN(ipt,l)                           &
     &                          / max(1.0e-10_kp,QLCN(ipt,l)+QICN(ipt,l))
              else
                QLCN(ipt,l)     = qli(l)
                QICN(ipt,l)     = qii(l)
                CNV_FICE(ipt,l) = qii(l)/max(1.0e-10_kp,qii(l)+qli(l))
              endif
!!            CNV_PRC3(ipt,l) = PCU(l)/dt
!             CNV_PRC3(ipt,l) = zero
!             if(PCU(l) < zero) write(*,*)"AAA777",PCU(l),ipt,l
              cf_upi(ipt,l)   = max(zero,min(0.02_kp*log(one+           &
     &                             500.0_kp*ud_mf(ipt,l)/dt), cfmax))
!    &                             500*ud_mf(ipt,l)/dt), 0.60))
              CLCN(ipt,l)     = cf_upi(ipt,l)  !downdraft is below updraft
              w_upi(ipt,l)    = ud_mf(ipt,l)*toi(l)*rgas /              &
     &                        (dt*max(cf_upi(ipt,l),1.0e-12_kp)*prsl(ipt,l))
            endif

            if (ntr > 0) then
              do n=1,ntr
                ccin(ipt,l,n+2) = uvi(l,n)           ! Tracers
              enddo
            endif
          enddo
          if (advcld) then
            do l=1,k
              ccin(ipt,l,1) = qii(l)          ! Cloud ice
              ccin(ipt,l,2) = qli(l)          ! Cloud water
            enddo
          else
            do l=1,k
              ccin(ipt,l,1) = ccin(ipt,l,1) + cli(l)
              ccin(ipt,l,2) = ccin(ipt,l,2) + clw(l)
            enddo
          endif
        endif
!
!     Velocity scale from the downdraft!
!
        DDVEL(ipt) = DDVEL(ipt) * DDFAC * GRAV / (prs(KP1)-prs(K))
!
      ENDDO                            ! End of the IPT Loop!

      deallocate (alfint, uvi, trcfac, rcu)
!
      RETURN
      end subroutine rascnv_run

!>\ingroup rascnv_schm
      SUBROUTINE CLOUD(                                                 &
     &                  K, KP1, KD, NTRC, KBLMX, kblmn                  &
     &,                 FRACBL, MAX_NEG_BOUY, vsmooth, aw_scal          &
     &,                 REVAP, WRKFUN, CALKBL, CRTFUN                   &
     &,                 DT, KDT, TLA, DPD                               &
     &,                 ALFINT, RHFACL, RHFACS, area, ccwf, cd, trcfac  &
     &,                 alfind, rhc_ls, phil, phih, prs, prsm, sgcs     &
     &,                 TOI, QOI, ROI,  QLI, QII, KPBL, DSFC            &
     &,                 TCU, QCU, RCU, PCU, FLX, FLXD, CUP, WFNC,fscav_ &
     &,                 trcmin, ntk, c0, qw0, c0i, qi0, dlq_fac)
!    &,                 ctei)

!
!***********************************************************************
!******************** Relaxed  Arakawa-Schubert ************************
!****************** Plug Compatible Scalar Version *********************
!************************ SUBROUTINE CLOUD  ****************************
!************************  October 2004     ****************************
!********************  VERSION 2.0  (modified) *************************
!************* Shrinivas.Moorthi@noaa.gov (301) 683-3718  ***** ********
!***********************************************************************
!*References:
!-----------
!     NOAA Technical Report NWS/NCEP 99-01:
!     Documentation of Version 2 of Relaxed-Arakawa-Schubert
!     Cumulus Parameterization with Convective Downdrafts, June 1999.
!     by S. Moorthi and M. J. Suarez.
!
!     Relaxed Arakawa-Schubert Cumulus Parameterization (Version 2)
!     with Convective Downdrafts - Unpublished Manuscript (2002)
!     by Shrinivas Moorthi and Max J. Suarez.
!
!***********************************************************************
!
!===>    UPDATES CLOUD TENDENCIES DUE TO A SINGLE CLOUD
!===>    DETRAINING AT LEVEL KD.
!
!***********************************************************************
!
!===>  TOI(K)     INOUT   TEMPERATURE             KELVIN
!===>  QOI(K)     INOUT   SPECIFIC HUMIDITY       NON-DIMENSIONAL
!===>  ROI(K,NTRC)INOUT   TRACER                  ARBITRARY
!===>  QLI(K)     INOUT   LIQUID WATER            NON-DIMENSIONAL
!===>  QII(K)     INOUT   ICE                     NON-DIMENSIONAL

!===>  PRS(KP1)   INPUT   PRESSURE @ EDGES        MB
!===>  PRSM(K)    INPUT   PRESSURE @ LAYERS       MB
!===>  SGCS(K)    INPUT   Local sigma
!===>  PHIH(KP1)  INPUT   GEOPOTENTIAL @ EDGES  IN MKS units
!===>  PHIL(K)    INPUT   GEOPOTENTIAL @ LAYERS IN MKS units
!===>  PRJ(KP1)   INPUT   (P/P0)^KAPPA  @ EDGES   NON-DIMENSIONAL
!===>  PRJM(K)    INPUT   (P/P0)^KAPPA  @ LAYERS  NON-DIMENSIONAL

!===>  K          INPUT   THE RISE & THE INDEX OF THE SUBCLOUD LAYER
!===>  KD         INPUT   DETRAINMENT LEVEL ( 1<= KD < K )          
!===>  NTRC       INPUT   NUMBER OF TRACERS. MAY BE ZERO.
!===>  kblmx      INPUT   highest level the pbl can take
!===>  kblmn      INPUT   lowest  level the pbl can take
!===>  DPD        INPUT   Critical normalized pressure (i.e. sigma) at the cloud top
!                         No downdraft calculation if the cloud top pressure is higher
!                         than DPD*PRS(KP1)
!
!===>  TCU(K  )   UPDATE  TEMPERATURE TENDENCY       DEG
!===>  QCU(K  )   UPDATE  WATER VAPOR TENDENCY       (G/G)
!===>  RCU(K,NTRC)UPDATE  TRACER TENDENCIES          ND
!===>  PCU(K)     UPDATE  PRECIP @ BASE OF LAYER     KG/M^2
!===>  FLX(K  )   UPDATE  MASS FLUX @ TOP OF LAYER   KG/M^2
!===>  CUP        UPDATE  PRECIPITATION AT THE SURFACE KG/M^2
!
      IMPLICIT NONE
!
      real (kind=kind_phys), parameter :: RHMAX=1.0_kp                   & ! MAX RELATIVE HUMIDITY
     &,                                   QUAD_LAM=1.0_kp                & ! MASK FOR QUADRATIC LAMBDA
     &,                                   RHRAM=0.05_kp                  & ! PBL RELATIVE HUMIDITY RAMP
!    &,                                   RHRAM=0.15_kp                 !& ! PBL RELATIVE HUMIDITY RAMP
     &,                                   HCRITD=4000.0_kp               & ! Critical Moist Static Energy for Deep clouds
!    &,                                   HCRITS=2000.0_kp               & ! Critical Moist Static Energy for Shallow clouds
     &,                                   HCRITS=2500.0_kp               & ! Critical Moist Static Energy for Shallow clouds
     &,                                   pcrit_lcl=250.0_kp             & ! Critical pressure difference between boundary layer top
                                                                          ! layer top and lifting condensation level (hPa)
!    &,                                   hpert_fac=1.01_kp             !& ! Perturbation on hbl when ctei=.true.
!    &,                                   hpert_fac=1.005_kp            !& ! Perturbation on hbl when ctei=.true.
     &,                                   qudfac=quad_lam*half           &
     &,                                   shalfac=3.0_kp                 &
!    &,                                   qudfac=quad_lam*pt25, shalfac=3.0_kp       !& !  Yogesh's
!    &,                                   c0ifac=0.07_kp                 & ! following Han et al, 2016 MWR
!    &,                                   c0ifac=0.001_kp                & ! following Han et al, 2017 Weather and Forecasting
     &,                                   c0ifac=0.0_kp                  &
     &,                                   dpnegcr  = 150.0_kp
!    &,                                   dpnegcr  = 100.0_kp
!    &,                                   dpnegcr  = 200.0_kp
!
      real(kind=kind_phys), parameter :: ERRMIN=0.0001_kp                &
     &,                                  ERRMI2=0.1_kp*ERRMIN            &
!    &,                                  rainmin=1.0e-9_kp              !&
     &,                                  rainmin=1.0e-8_kp               &
     &,                                  oneopt9=one/0.09_kp             &
     &,                                  oneopt4=one/0.04_kp
      real(kind=kind_phys), parameter :: almax=1.0e-2_kp                 &
     &,                                  almin1=0.0_kp,   almin2=0.0_kp
      real(kind=kind_phys), parameter :: bldmax=300.0_kp, bldmin=25.0_kp
!
!  INPUT ARGUMENTS

!     LOGICAL REVAP, WRKFUN, CALKBL, CRTFUN, CALCUP, ctei
      LOGICAL REVAP, WRKFUN, CALKBL, CRTFUN, CALCUP
      logical vsmooth, aw_scal
      INTEGER K, KP1, KD, NTRC, kblmx, kblmn, ntk


      real(kind=kind_phys), dimension(K)   ::  TOI,  QOI, PRSM, QLI, QII&
     &,                                        PHIL, SGCS, rhc_ls       &
     &,                                        alfind
      real(kind=kind_phys), dimension(KP1) ::  PRS,  PHIH
      real(kind=kind_phys), dimension(K,NTRC) :: ROI, trcfac
      real(kind=kind_phys), dimension(ntrc)   :: trcmin
      real(kind=kind_phys)                    :: CD, DSFC
      INTEGER                                 :: KPBL, KBL, KB1, kdt

      real(kind=kind_phys) ALFINT(K,NTRC+4)
      real(kind=kind_phys) FRACBL, MAX_NEG_BOUY, DPD                    &
     &,                    RHFACL, RHFACS, area, ccwf                   &
     &,                    c0, qw0, c0i, qi0, dlq_fac
 
!  UPDATE ARGUMENTS

      real(kind=kind_phys), dimension(K)      :: TCU, QCU, TCD, QCD, PCU
      real(kind=kind_phys), dimension(KP1)    :: FLX, FLXD
      real(kind=kind_phys), dimension(K,NTRC) :: RCU
      real(kind=kind_phys)                    :: CUP
!
!  TEMPORARY WORK SPACE

      real(kind=kind_phys), dimension(KD:K)   :: HOL, QOL, HST, QST     &
     &,                           TOL, GMH, AKT, AKC, BKC, LTL, RNN     &
     &,                           FCO, PRI, QIL, QLL, ZET, XI, RNS      &
     &,                           Q0U, Q0D, vtf, CIL, CLL, ETAI, dlq    &
     &,                           wrk1, wrk2, dhdp, qrb, qrt, evp       &
     &,                           ghd, gsd, etz, cldfr,  sigf, rho

      real(kind=kind_phys), dimension(KD:KP1) :: GAF, GMS, GAM, DLB     &
     &,                           DLT, ETA, PRL, BUY, ETD, HOD, QOD, wvl
      real(kind=kind_phys), dimension(KD:K-1) :: etzi

      real(kind=kind_phys) fscav_(ntrc)

      LOGICAL ep_wfn, cnvflg, LOWEST, DDFT, UPDRET

      real(kind=kind_phys) ALM,   DET,    HCC,  CLP                     &
     &,                    HSU,   HSD,    QTL,  QTV                     &
     &,                    AKM,   WFN,    HOS,  QOS                     &
     &,                    AMB,   TX1,    TX2,  TX3                     &
     &,                    TX4,   TX5,    QIS,  QLS                     &
     &,                    HBL,   QBL,    RBL(NTRC), wcbase             &
     &,                    QLB,   QIB,    PRIS                          &
     &,                    WFNC,  TX6,    ACR                           &
     &,                    TX7,   TX8,    TX9,  RHC                     &
     &,                    hstkd, qstkd,  ltlkd, q0ukd, q0dkd, dlbkd    &
     &,                    qtp, qw00, qi00, qrbkd                       &
     &,                    hstold, rel_fac, prism                       &
     &,                    TL, PL, QL, QS, DQS, ST1, SGN, TAU,          &
     &                     QTVP, HB, QB, TB, QQQ,                       &
     &                     HCCP, DS, DH, AMBMAX, X00, EPP, QTLP,        &
     &                     DPI, DPHIB, DPHIT, DEL_ETA, DETP,            &
     &                     TEM, TEM1, TEM2, TEM3, TEM4,                 &
     &                     ST2, ST3, ST4, ST5,                          &
     &                     ERRH, ERRW, ERRE, TEM5,                      &
     &                     TEM6, HBD, QBD, st1s, shal_fac, hmax, hmin,  &
     &                     dhdpmn, avt, avq, avr, avh                   &
     &,                    TRAIN, DOF, CLDFRD, tla, gmf                 &
     &,                    FAC, RSUM1, RSUM2, RSUM3, dpneg, hcrit       &
     &,                    ACTEVAP,AREARAT,DELTAQ,MASS,MASSINV,POTEVAP  &
     &,                    TEQ,QSTEQ,DQDT,QEQ                           &
     &,                    CLFRAC, DT,      clvfr, delzkm, fnoscav, delp
!    &,                    CLFRAC, DT, clf, clvfr, delzkm, fnoscav, delp
!    &,                    almin1, almin2

      INTEGER I, L,  N,  KD1, II, iwk, idh, lcon                        &
     &,       IT, KM1, KTEM, KK, KK1, LM1, LL, LP1, kbls, kmxh          &
     &,       kblh, kblm, kblpmn, kmax, kmaxm1, kmaxp1, klcl, kmin, kmxb
!
!***********************************************************************
!
!     almin2 = 0.2 * sqrt(pi/area)
!     almin1 = almin2

      KM1 = K  - 1
      KD1 = KD + 1

      do l=1,K
        tcd(L) = zero
        qcd(L) = zero
      enddo
!
      CLDFRD   = zero
      DOF      = zero
      PRL(KP1) = PRS(KP1)
!
      DO L=KD,K
        RNN(L)   = zero
        ZET(L)   = zero
        XI(L)    = zero
!
        TOL(L)   = TOI(L)
        QOL(L)   = QOI(L)
        PRL(L)   = PRS(L)
        CLL(L)   = QLI(L)
        CIL(L)   = QII(L)
        BUY(L)   = zero

        wvl(l)   = zero
      ENDDO
      wvl(kp1) = zero
!
      if (vsmooth) then
        do l=kd,k
          wrk1(l) = tol(l)
          wrk2(l) = qol(l)
        enddo
        do l=kd1,km1
          tol(l) = pt25*wrk1(l-1) + half*wrk1(l) + pt25*wrk1(l+1)
          qol(l) = pt25*wrk2(l-1) + half*wrk2(l) + pt25*wrk2(l+1)
        enddo
      endif
!
      DO L=KD, K
        DPI    = ONE / (PRL(L+1) - PRL(L))
        PRI(L) = GRAVFAC * DPI
!
        PL     = PRSM(L)
        TL     = TOL(L)

        rho(l) = cmb2pa * pl / (rgas*tl*(one+nu*qol(l)))

        AKT(L) = (PRL(L+1) - PL) * DPI
!
        CALL QSATCN(TL, PL, QS, DQS)
!
        QST(L) = QS
        GAM(L) = DQS * ELOCP
        ST1    = ONE + GAM(L)
        GAF(L) = ONEOALHL * GAM(L) / ST1
 
        QL     = MAX(MIN(QS*RHMAX,QOL(L)), ONE_M10)
        QOL(L) = QL
 
        TEM    = CP * TL
        LTL(L) = TEM * ST1 / (ONE+NU*(QST(L)+TL*DQS))
        vtf(L) = one + NU * QL
        ETA(L) = ONE / (LTL(L) * VTF(L))

        HOL(L) = TEM + QL * ALHL
        HST(L) = TEM + QS * ALHL
!
      ENDDO
!
      ETA(KP1) = ZERO
      GMS(K)   = ZERO
!
      AKT(KD)  = HALF
      GMS(KD)  = ZERO
!
      CLP      = ZERO
!
      GAM(KP1) = GAM(K)
      GAF(KP1) = GAF(K)
!
      DO L=K,KD1,-1
        DPHIB  = PHIL(L) - PHIH(L+1)
        DPHIT  = PHIH(L) - PHIL(L)
!
        DLB(L) = DPHIB * ETA(L) ! here eta contains 1/(L*(1+nu*q))
        DLT(L) = DPHIT * ETA(L)
!
        QRB(L) = DPHIB
        QRT(L) = DPHIT
!
        ETA(L) = ETA(L+1) + DPHIB

        HOL(L) = HOL(L) + ETA(L)
        hstold = hst(l)
        HST(L) = HST(L) + ETA(L)
!
        ETA(L) = ETA(L) + DPHIT
      ENDDO
!
!     For the cloud top layer
!
      L = KD

      DPHIB  = PHIL(L) - PHIH(L+1)
!
      DLB(L) = DPHIB * ETA(L)
!
      QRB(L) = DPHIB
      QRT(L) = DPHIB
!
      ETA(L) = ETA(L+1) + DPHIB

      HOL(L) = HOL(L) + ETA(L)
      HST(L) = HST(L) + ETA(L)
!
!     To determine KBL internally -- If KBL is defined externally
!     the following two loop should be skipped
!
      if (sgcs(kd) < 0.5_kp) then
         hcrit = hcritd
      elseif (sgcs(kd) > 0.65_kp) then
         hcrit = hcrits
      else
         hcrit = (hcrits*(sgcs(kd)-0.5_kp) + hcritd*(0.65_kp-sgcs(kd)))&
     &         * (one/0.15_kp)
      endif
      IF (CALKBL) THEN
         KTEM = MAX(KD+1, KBLMX)
         hmin = hol(k)
         kmin = k
         do l=km1,kd,-1
           if (hmin > hol(l)) then
             hmin = hol(l)
             kmin = l
           endif
         enddo
         if (kmin == k) return
         hmax = hol(k)
         kmax = k
         do l=km1,ktem,-1
           if (hmax < hol(l)) then
             hmax = hol(l)
             kmax = l
           endif
         enddo
         kmxb = kmax
         if (kmax < kmin) then
           kmax = k
           kmxb = k
           hmax = hol(kmax)
         elseif (kmax < k) then
           do l=kmax+1,k
!            if (abs(hol(kmax)-hol(l)) > half * hcrit) then
             if (abs(hol(kmax)-hol(l)) >        hcrit) then
               kmxb = l - 1
               exit
             endif
           enddo
         endif
         kmaxm1 = kmax - 1
         kmaxp1 = kmax + 1
         kblpmn = kmax
!
         dhdp(kmax:k) = zero
         dhdpmn = dhdp(kmax)
         do l=kmaxm1,ktem,-1
           dhdp(l) = (HOL(L)-HOL(L+1)) / (PRL(L+2)-PRL(L))
           if (dhdp(l) < dhdpmn) then
             dhdpmn = dhdp(l)
             kblpmn = l + 1
           elseif (dhdp(l) > zero .and. l <= kmin) then
             exit
           endif
         enddo
         kbl = kmax
         if (kblpmn < kmax) then
           do l=kblpmn,kmaxm1
             if (hmax-hol(l) < half*hcrit) then
               kbl = l
               exit
             endif
           enddo
         endif
!
         klcl = kd1
         if (kmax > kd1) then
           do l=kmaxm1,kd1,-1
             if (hmax > hst(l)) then
               klcl = l+1
               exit
             endif
           enddo
         endif

!        if (klcl == kd .or. klcl < ktem) return

!        This is to handle mid-level convection from quasi-uniform h

         if (kmax < kmxb) then
           kmax   = max(kd1, min(kmxb,k))
           kmaxm1 = kmax - 1
           kmaxp1 = kmax + 1
         endif


!        if (prl(Kmaxp1) - prl(klcl) > 250.0 ) return

         ii  = max(kbl,kd1)
         kbl = max(klcl,kd1)
         tem = min(50.0_kp,max(10.0_kp,(prl(kmaxp1)-prl(kd))*0.10_kp))
         if (prl(kmaxp1) - prl(ii) > tem .and. ii > kbl) kbl = ii

         if (kbl .ne. ii) then
           if (PRL(kmaxp1)-PRL(KBL) > bldmax) kbl = max(kbl,ii)
         endif
         if (kbl < ii) then
           if (hol(ii)-hol(ii-1) > half*hcrit) kbl = ii
         endif

         if (prl(kbl) - prl(klcl) > pcrit_lcl) return
!
!        KBL  = min(kmax, MAX(KBL,KBLMX))
         KBL  = min(kblmn, MAX(KBL,KBLMX))
!        kbl  = min(kblh,kbl)
!!!
!        tem1 = max(prl(kP1)-prl(k),                                    &
!    &                     min((prl(kbl) - prl(kd))*0.05, 10.0))
!!   &                     min((prl(kbl) - prl(kd))*0.05, 20.0))
!!   &                     min((prl(kbl) - prl(kd))*0.05, 30.0))
!        if (prl(kp1)-prl(kbl) < tem1) then
!          KTEM = MAX(KD+1, KBLMX)
!          do l=k,KTEM,-1
!            tem = prl(kp1) - prl(l)
!            if (tem > tem1) then
!              kbl = min(kbl,l)
!              exit
!            endif
!          enddo
!        endif
!        if (kbl == kblmx .and. kmax >= km1) kbl = k - 1
!!!

         KPBL = KBL

      ELSE
         KBL = KPBL
      ENDIF
!
      KBL = min(kmax,MAX(KBL,KD+2))
      KB1 = KBL - 1
!
      if (PRL(Kmaxp1)-PRL(KBL) > bldmax .or. kb1 <= kd ) then
!    &          .or. PRL(Kmaxp1)-PRL(KBL) < bldmin) then
        return
      endif
!
      PRIS     = ONE / (PRL(KP1)-PRL(KBL))
      PRISM    = ONE / (PRL(Kmaxp1)-PRL(KBL))
      TX1      = ETA(KBL)         ! geopotential height at KBL
!
      GMS(KBL) = zero
      XI(KBL)  = zero
      ZET(KBL) = zero
!
      shal_fac = one
!     if (prl(kbl)-prl(kd) < 300.0 .and. kmax == k) shal_fac = shalfac
      if (prl(kbl)-prl(kd) < 350.0_kp .and. kmax == k) shal_fac = shalfac
      DO L=Kmax,KD,-1
        IF (L >= KBL) THEN
          ETA(L) = (PRL(Kmaxp1)-PRL(L)) * PRISM
        ELSE
          ZET(L) = (ETA(L) - TX1) * ONEBG
          XI(L)  =  ZET(L) * ZET(L) * (QUDFAC*shal_fac)
          ETA(L) =  ZET(L) - ZET(L+1)
          GMS(L) =  XI(L)  - XI(L+1)
        ENDIF
      ENDDO
      if (kmax < k) then
        do l=kmaxp1,kp1
          eta(l) = zero
        enddo
      endif
!
      HBL = HOL(Kmax) * ETA(Kmax)
      QBL = QOL(Kmax) * ETA(Kmax)
      QLB = CLL(Kmax) * ETA(Kmax)
      QIB = CIL(Kmax) * ETA(Kmax)
      TX1 = QST(Kmax) * ETA(Kmax)
!
      DO L=Kmaxm1,KBL,-1
         TEM = ETA(L) - ETA(L+1)
         HBL = HBL + HOL(L) * TEM
         QBL = QBL + QOL(L) * TEM
         QLB = QLB + CLL(L) * TEM
         QIB = QIB + CIL(L) * TEM
         TX1 = TX1 + QST(L) * TEM
      ENDDO

!     if (ctei .and. sgcs(kd) > 0.65) then
!        hbl = hbl * hpert_fac
!        qbl = qbl * hpert_fac
!     endif

!                                   Find Min value of HOL in TX2
      TX2 = HOL(KD)
      IDH = KD1
      DO L=KD1,KB1
        IF (HOL(L) < TX2) THEN
           TX2 = HOL(L)
           IDH = L             ! Level of minimum moist static energy!
        ENDIF
      ENDDO
      IDH = 1
!     IDH = MAX(KD1, IDH)
      IDH = MAX(KD, IDH)       ! Moorthi May, 31, 2019
!
      TEM1 = HBL - HOL(KD)
      TEM  = HBL - HST(KD1) - LTL(KD1) * NU *(QOL(KD1)-QST(KD1))
      LOWEST = KD == KB1

      lcon = kd
      do l=kb1,kd1,-1
        if (hbl >= hst(l)) then
          lcon = l
          exit
        endif
      enddo
!
      if (lcon == kd .or. kbl <= kd .or. prl(kbl)-prsm(lcon) > 150.0_kp) &
     &                                    return
!
      TX1    = RHFACS - QBL / TX1       !     Average RH

      cnvflg = (TEM > ZERO .OR. (LOWEST .AND. TEM1 >= ZERO))            &
     &         .AND. TX1 < RHRAM

      IF (.NOT. cnvflg) RETURN
!
      RHC = MAX(ZERO, MIN(ONE, EXP(-20.0_kp*TX1) ))
!
      wcbase = 0.1_kp
      if (ntrc > 0) then
        DO N=1,NTRC
          RBL(N) = ROI(Kmax,N) * ETA(Kmax)
        ENDDO
        DO N=1,NTRC
          DO L=KmaxM1,KBL,-1
            RBL(N) = RBL(N) + ROI(L,N)*(ETA(L)-ETA(L+1))
          ENDDO
        ENDDO
!
!       if (ntk > 0 .and. aw_scal) then
        if (ntk > 0) then
          if (rbl(ntk) > zero) then
            wcbase = min(two, max(wcbase, sqrt(twoo3*rbl(ntk))))
!           wcbase = min(one, max(wcbase, sqrt(twoo3*rbl(ntk))))
          endif
        endif

      endif
!
      TX4 = zero
      TX5 = zero
!
      TX3 = QST(KBL) - GAF(KBL) * HST(KBL)
      DO L=KBL,K
        QIL(L) = MAX(ZERO, MIN(ONE, (TCR-TCL-TOL(L))*TCRF))
      ENDDO
!
      DO L=KB1,KD1,-1
        lp1      = l + 1
        TEM      = QST(L) - GAF(L) * HST(L)
        TEM1     = (TX3 + TEM) * half
        ST2      = (GAF(L)+GAF(LP1)) * half
!
        FCO(LP1) =            TEM1 + ST2 * HBL

        RNN(LP1) = ZET(LP1) * TEM1 + ST2 * TX4
        GMH(LP1) = XI(LP1)  * TEM1 + ST2 * TX5
!
        TX3      = TEM
        TX4      = TX4 + ETA(L) * HOL(L)
        TX5      = TX5 + GMS(L) * HOL(L)
!
        QIL(L)   = MAX(ZERO, MIN(ONE, (TCR-TCL-TOL(L))*TCRF))
        QLL(LP1) = (half*ALHF) * ST2 * (QIL(L)+QIL(LP1)) + ONE
      ENDDO
!
!     FOR THE CLOUD TOP -- L=KD
!
      L = KD
!
      lp1      = l + 1
      TEM      = QST(L) - GAF(L) * HST(L)
      TEM1     = (TX3 + TEM) * half
      ST2      = (GAF(L)+GAF(LP1)) * half
!
      FCO(LP1) =            TEM1 + ST2 * HBL
      RNN(LP1) = ZET(LP1) * TEM1 + ST2 * TX4
      GMH(LP1) = XI(LP1)  * TEM1 + ST2 * TX5
!
      FCO(L)   = TEM + GAF(L) * HBL
      RNN(L)   = TEM * ZET(L) + (TX4 + ETA(L)*HOL(L)) * GAF(L)
      GMH(L)   = TEM * XI(L)  + (TX5 + GMS(L)*HOL(L)) * GAF(L)
!
!   Replace FCO for the Bottom
!
      FCO(KBL) = QBL
      RNN(KBL) = zero
      GMH(KBL) = zero
!
      QIL(KD)  =  MAX(ZERO, MIN(ONE, (TCR-TCL-TOL(KD))*TCRF))
      QLL(KD1) = (half*ALHF) * ST2 * (QIL(KD) + QIL(KD1)) + ONE
      QLL(KD ) = ALHF * GAF(KD) * QIL(KD) + ONE
!
      st1  = qil(kd)
      st2  = c0i * st1
      if (c0ifac > 1.0e-6_kp) st2  = st2 * exp(c0ifac*min(tol(kd)-t0c,zero))
      tem  = c0  * (one-st1)
      tem2 = st2*qi0 + tem*qw0
!
      DO L=KD,KB1
         lp1    = l + 1
         tx2    = akt(l) * eta(l)
         tx1    = tx2 * tem2
         q0u(l) = tx1
         FCO(L) = FCO(LP1) - FCO(L) + tx1
         RNN(L) = RNN(LP1) - RNN(L)                                     &
     &          + ETA(L)*(QOL(L)+CLL(L)+CIL(L)) + tx1*zet(l)
         GMH(L) = GMH(LP1) - GMH(L)                                     &
     &          + GMS(L)*(QOL(L)+CLL(L)+CIL(L)) + tx1*xi(l)
!
         tem1   = (one-akt(l)) * eta(l)

         AKT(L) = QLL(L)   + (st2 + tem) * tx2

         AKC(L) = one / AKT(L)
!
         st1    = half * (qil(l)+qil(lp1))
         st2    = c0i * st1
         if (c0ifac > 1.0e-6_kp) st2 = st2 * exp(c0ifac*min(tol(lp1)-t0c,zero))
         tem    = c0  * (one-st1)
         tem2   = st2*qi0 + tem*qw0
!
         BKC(L) = QLL(LP1) - (st2 + tem) * tem1
!
         tx1    = tem1*tem2
         q0d(l) = tx1
         FCO(L) = FCO(L) + tx1
         RNN(L) = RNN(L) + tx1*zet(lp1)
         GMH(L) = GMH(L) + tx1*xi(lp1)
      ENDDO

      qw00 = qw0
      qi00 = qi0
      ii = 0
  777 continue
!
      ep_wfn = .false.
      RNN(KBL) = zero
      TX3      = bkc(kb1) * (QIB + QLB)
      TX4      = zero
      TX5      = zero
      DO L=KB1,KD1,-1
        TEM    = BKC(L-1)       * AKC(L)
        TX3    = (TX3 + FCO(L)) * TEM
        TX4    = (TX4 + RNN(L)) * TEM
        TX5    = (TX5 + GMH(L)) * TEM
      ENDDO
      IF (KD < KB1) THEN
         HSD   = HST(KD1) + LTL(KD1) *  NU *(QOL(KD1)-QST(KD1))
      ELSE
         HSD   = HBL
      ENDIF
!
      TX3 = (TX3 + FCO(KD)) * AKC(KD)
      TX4 = (TX4 + RNN(KD)) * AKC(KD)
      TX5 = (TX5 + GMH(KD)) * AKC(KD)
      ALM = ALHF*QIL(KD) - LTL(KD) * VTF(KD)
!
      HSU = HST(KD) + LTL(KD) * NU * (QOL(KD)-QST(KD))
!
!===> VERTICAL INTEGRALS NEEDED TO COMPUTE THE ENTRAINMENT PARAMETER
!
      TX1 = ALM * TX4
      TX2 = ALM * TX5

      DO L=KD,KB1
        TAU = HOL(L) - HSU
        TX1 = TX1 + TAU * ETA(L)
        TX2 = TX2 + TAU * GMS(L)
      ENDDO
!
!     MODIFY HSU TO INCLUDE CLOUD LIQUID WATER AND ICE TERMS
!
      HSU    = HSU - ALM * TX3
!
      CLP    = ZERO
      ALM    = -100.0_kp
      HOS    = HOL(KD)
      QOS    = QOL(KD)
      QIS    = CIL(KD)
      QLS    = CLL(KD)

      cnvflg = HBL > HSU .and. abs(tx1) > 1.0e-4_kp

!***********************************************************************

      ST1 = HALF*(HSU + HSD)

      IF (cnvflg) THEN
!
!  STANDARD CASE:
!   CLOUD CAN BE NEUTRALLY BOUYANT AT MIDDLE OF LEVEL KD W/ +VE LAMBDA.
!   EPP < .25 IS REQUIRED TO HAVE REAL ROOTS.
!
        clp = one
        st2 = hbl - hsu
!
        if (tx2 == zero) then
          alm = - st2 / tx1
          if (alm > almax) alm = -100.0_kp
        else
          x00 = tx2 + tx2
          epp = tx1 * tx1 - (x00+x00)*st2
          if (epp > zero) then
            x00  = one / x00
            tem  = sqrt(epp)
            tem1 = (-tx1-tem)*x00
            tem2 = (-tx1+tem)*x00
            if (tem1 > almax) tem1 = -100.0_kp
            if (tem2 > almax) tem2 = -100.0_kp
            alm  = max(tem1,tem2)

          endif
        endif
!
!  CLIP CASE:
!   NON-ENTRAINIG CLOUD DETRAINS IN LOWER HALF OF TOP LAYER.
!   NO CLOUDS ARE ALLOWED TO DETRAIN BELOW THE TOP LAYER.
!
      ELSEIF (HBL <= HSU .AND. HBL > ST1) THEN
        ALM = ZERO
!       CLP = (HBL-ST1) / (HSU-ST1)    ! commented on Jan 16, 2010
      ENDIF
!
      cnvflg = .TRUE.
      IF (ALMIN1 > zero) THEN
        IF (ALM >= ALMIN1) cnvflg = .FALSE.
      ELSE
        LOWEST = KD == KB1
        IF ( (ALM > ZERO) .OR.                                          &
     &     (.NOT. LOWEST .AND. ALM == ZERO) ) cnvflg = .FALSE.
      ENDIF
!
!===>  IF NO SOUNDING MEETS SECOND CONDITION, RETURN
!
      IF (cnvflg) THEN
         IF (ii > 0 .or. (qw00 == zero .and. qi00 == zero)) RETURN
         CLP = one
         ep_wfn = .true.
         GO TO 888
      ENDIF
!
      st1s = ONE
      IF(CLP > ZERO .AND. CLP < ONE) THEN
        ST1     = HALF*(ONE+CLP)
        ST2     = ONE - ST1
        st1s    = st1
        hstkd   = hst(kd)
        qstkd   = qst(kd)
        ltlkd   = ltl(kd)
        q0ukd   = q0u(kd)
        q0dkd   = q0d(kd)
        dlbkd   = dlb(kd)
        qrbkd   = qrb(kd)
!
        HST(KD) = HST(KD)*ST1 + HST(KD1)*ST2
        HOS     = HOL(KD)*ST1 + HOL(KD1)*ST2
        QST(KD) = QST(KD)*ST1 + QST(KD1)*ST2
        QOS     = QOL(KD)*ST1 + QOL(KD1)*ST2
        QLS     = CLL(KD)*ST1 + CLL(KD1)*ST2
        QIS     = CIL(KD)*ST1 + CIL(KD1)*ST2
        LTL(KD) = LTL(KD)*ST1 + LTL(KD1)*ST2
!
        DLB(KD) = DLB(KD)*CLP
        qrb(KD) = qrb(KD)*CLP
        ETA(KD) = ETA(KD)*CLP
        GMS(KD) = GMS(KD)*CLP
        Q0U(KD) = Q0U(KD)*CLP
        Q0D(KD) = Q0D(KD)*CLP
      ENDIF
!
!
!***********************************************************************
!
!    Critical workfunction is included in this version
!
      ACR  = zero
      TEM  = PRL(KD1) - (PRL(KD1)-PRL(KD)) * CLP * HALF
      tx1  = PRL(KBL) - TEM
      tx2  = min(900.0_kp, max(tx1,100.0_kp))
      tem1 = log(tx2*0.01_kp) * oneolog10
      tem2 = one - tem1
      if ( kdt == 1 ) then
!       rel_fac = (dt * facdt)  / (tem1*12.0_kp + tem2*3.0)
        rel_fac = (dt * facdt)  / (tem1*6.0_kp + tem2*adjts_s)
      else
        rel_fac = (dt * facdt) / (tem1*adjts_d + tem2*adjts_s)
      endif
!
!     rel_fac = max(zero, min(one,rel_fac))
      rel_fac = max(zero, min(half,rel_fac))
      
      IF (CRTFUN) THEN
        iwk = tem*0.02_kp - 0.999999999_kp
        iwk = MAX(1, MIN(iwk, 16))
        ACR = tx1 * (AC(iwk) + tem * AD(iwk)) * CCWF
      ENDIF
!
!===>  NORMALIZED MASSFLUX
!
!  ETA IS THE THICKNESS COMING IN AND normalized  MASS FLUX GOING OUT.
!  GMS IS THE THICKNESS SQUARE ; IT IS LATER REUSED FOR GAMMA_S
!
!     ETA(K) = ONE

      DO L=KB1,KD,-1
        ETA(L)  = ETA(L+1) + ALM * (ETA(L) + ALM * GMS(L))
        ETAI(L) = one / ETA(L)
      ENDDO
      ETAI(KBL) = one
!
!===>  CLOUD WORKFUNCTION
!
      WFN    = ZERO
      AKM    = ZERO
      DET    = ZERO
      HCC    = HBL
      cnvflg = .FALSE.
      QTL    = QST(KB1) - GAF(KB1)*HST(KB1)
      TX1    = HBL
!
      qtv    = qbl
      det    = qlb + qib
!
      tx2    = zero
      dpneg  = zero
!
      DO L=KB1,KD1,-1
         lm1 = l - 1
         lp1 = l + 1
         DEL_ETA = ETA(L) - ETA(LP1)
         HCCP = HCC + DEL_ETA*HOL(L)
!
         QTLP = QST(LM1) - GAF(LM1)*HST(LM1)
         QTVP = half * ((QTLP+QTL)*ETA(L)                               &
     &               + (GAF(L)+GAF(LM1))*HCCP)
         ST1  = ETA(L)*Q0U(L) + ETA(LP1)*Q0D(L)
         DETP = (BKC(L)*DET - (QTVP-QTV)                                &
     &        + DEL_ETA*(QOL(L)+CLL(L)+CIL(L)) + ST1)  * AKC(L)

         TEM1   = AKT(L)   - QLL(L)
         TEM2   = QLL(LP1) - BKC(L)
         RNS(L) = TEM1*DETP  + TEM2*DET - ST1

         qtp    = half * (qil(L)+qil(LM1))
         tem2   = min(qtp*(detp-eta(l)*qw00),                           &
     &               (one-qtp)*(detp-eta(l)*qi00))
         st1    = min(tx2,tem2)
         tx2    = tem2
!
         IF (rns(l) < zero .or. st1 < zero) ep_wfn = .TRUE.
         IF (DETP <= ZERO) cnvflg = .TRUE.

         ST1  = HST(L) - LTL(L)*NU*(QST(L)-QOL(L))


         TEM2 = HCCP   + DETP   * QTP * ALHF
!
         ST2  = LTL(L) * VTF(L)
         TEM5 = CLL(L) + CIL(L)
         TEM3 = (TX1  - ETA(LP1)*ST1 - ST2*(DET-TEM5*eta(lp1))) * DLB(L)
         TEM4 = (TEM2 - ETA(L  )*ST1 - ST2*(DETP-TEM5*eta(l)))  * DLT(L)
!
         ST1  = TEM3 + TEM4

         WFN = WFN + ST1       
         AKM = AKM - min(ST1,ZERO)

         if (st1 < zero .and. wfn < zero) then
           dpneg = dpneg + prl(lp1) - prl(l)
         endif

         BUY(L) = half * (tem3/(eta(lp1)*qrb(l)) + tem4/(eta(l)*qrt(l)))
!
         HCC = HCCP
         DET = DETP
         QTL = QTLP
         QTV = QTVP
         TX1 = TEM2

      ENDDO

      DEL_ETA = ETA(KD) - ETA(KD1)
      HCCP    = HCC + DEL_ETA*HOS
!
      QTLP    = QST(KD) - GAF(KD)*HST(KD)
      QTVP    = QTLP*ETA(KD) + GAF(KD)*HCCP
      ST1     = ETA(KD)*Q0U(KD) + ETA(KD1)*Q0D(KD)
      DETP    = (BKC(KD)*DET - (QTVP-QTV)                               &
     &        + DEL_ETA*(QOS+QLS+QIS) + ST1) * AKC(KD)
!
      TEM1    = AKT(KD)  - QLL(KD)
      TEM2    = QLL(KD1) - BKC(KD)
      RNS(KD) = TEM1*DETP  + TEM2*DET - ST1
!
      IF (rns(kd) < zero) ep_wfn = .TRUE.
      IF (DETP   <= ZERO) cnvflg = .TRUE.
!
  888 continue

      if (ep_wfn) then
        IF ((qw00 == zero .and. qi00 == zero)) RETURN
        if (ii == 0) then
          ii  = 1
          if (clp > zero .and. clp < one) then
            hst(kd) = hstkd
            qst(kd) = qstkd
            ltl(kd) = ltlkd
            q0u(kd) = q0ukd
            q0d(kd) = q0dkd
            dlb(kd) = dlbkd
            qrb(kd) = qrbkd
          endif
          do l=kd,kb1
            lp1    = l + 1
            FCO(L) = FCO(L) - q0u(l) - q0d(l)
            RNN(L) = RNN(L) - q0u(l)*zet(l) - q0d(l)*zet(lp1)
            GMH(L) = GMH(L) - q0u(l)*xi(l)  - q0d(l)*zet(lp1)
            ETA(L) = ZET(L) - ZET(LP1)
            GMS(L) = XI(L)  - XI(LP1)
            Q0U(L) = zero
            Q0D(L) = zero
          ENDDO
          qw00 = zero
          qi00 = zero

          go to 777
        else
          cnvflg = .true.
        endif
      endif
!
!
!     ST1 = 0.5 * (HST(KD)  - LTL(KD)*NU*(QST(KD)-QOS)                  &
!    &          +  HST(KD1) - LTL(KD1)*NU*(QST(KD1)-QOL(KD1)))
!
      ST1  = HST(KD) - LTL(KD)*NU*(QST(KD)-QOS)
      ST2  = LTL(KD) * VTF(KD)
      TEM5 = (QLS + QIS) * eta(kd1)
      ST1  = HALF * (TX1-ETA(KD1)*ST1-ST2*(DET-TEM5))*DLB(KD)
!
      WFN  = WFN + ST1
      AKM  = AKM - min(ST1,ZERO)   ! Commented on 08/26/02 - does not include top
!

      BUY(KD) = ST1 / (ETA(KD1)*qrb(kd))
!
      DET = DETP
      HCC = HCCP
      AKM = AKM / WFN


!***********************************************************************
!
      IF (WRKFUN) THEN ! If only to calculate workfunction save it and return
        IF (WFN >= zero) WFNC = WFN
        RETURN
      ELSEIF (.NOT. CRTFUN) THEN
        ACR = WFNC
      ENDIF
!
!===>  THIRD CHECK BASED ON CLOUD WORKFUNCTION
!
      CALCUP = .FALSE.

      TEM = max(0.05_kp, MIN(CD*200.0_kp, MAX_NEG_BOUY))
      IF (.not. cnvflg    .and. WFN > ACR .and.                         &
     &    dpneg < dpnegcr .and. AKM <= TEM)     CALCUP = .TRUE.
!
!===>  IF NO SOUNDING MEETS THIRD CONDITION, RETURN
!
      IF (.NOT. CALCUP) RETURN
!
! This is for not LL - 20050601
!     IF (ALMIN2 .NE. zero) THEN
!       IF (ALMIN1 .NE. ALMIN2) ST1 = one / max(ONE_M10,(ALMIN2-ALMIN1))
!       IF (ALM < ALMIN2) THEN
!          CLP = CLP * max(zero, min(one,(0.3 + 0.7*(ALM-ALMIN1)*ST1)))
!!         CLP = CLP * max(0.0, min(1.0,(0.2 + 0.8*(ALM-ALMIN1)*ST1)))
!!         CLP = CLP * max(0.0, min(1.0,(0.1 + 0.9*(ALM-ALMIN1)*ST1)))
!       ENDIF
!     ENDIF
!
      CLP = CLP * RHC
      dlq = zero
      tem = one / (one + dlq_fac)
      do l=kd,kb1
        rnn(l) = rns(l) * tem
        dlq(l) = rns(l) * tem * dlq_fac
      enddo
      DO L=KBL,K 
        RNN(L) = zero
      ENDDO
!
!     If downdraft is to be invoked, do preliminary check to see
!     if enough rain is available and then call DDRFT.
!
      DDFT = .FALSE.
      IF (dpd > zero) THEN
        TRAIN = zero
        IF (CLP > zero) THEN
          DO L=KD,KB1
            TRAIN = TRAIN + RNN(L)
          ENDDO
        ENDIF

        PL = (PRL(KD1) + PRL(KD))*HALF
        IF (TRAIN > 1.0e-4_kp .AND. PL <= dpd*prl(kp1)) DDFT  = .TRUE.
      ENDIF
!
      IF (DDFT) THEN ! Downdraft scheme based on (Cheng and Arakawa, 1997)
        CALL DDRFT(                                                     &
     &              K,   KP1,    KD                                     &
     &,             TLA, ALFIND, wcbase                                 &
     &,             TOL, QOL, HOL,   PRL, QST, HST, GAM, GAF            &
!    &,             TOL, QOL, HOL,   PRL, QST, HST, GAM, GAF, HBL, QBL  &
     &,             QRB, QRT, BUY,   KBL, IDH, ETA, RNN, ETAI           &
     &,             ALM, WFN, TRAIN, DDFT                               &
     &,             ETD, HOD, QOD,   EVP, DOF, CLDFR, ETZ               &
     &,             GMS, GSD, GHD,   wvl)               

      ENDIF
!
!  No Downdraft case (including case with no downdraft solution)
!  ---------------------------------------------------------
!
      IF (.NOT. DDFT) THEN
        DO L=KD,KP1
          ETD(L) = zero
          HOD(L) = zero
          QOD(L) = zero
          wvl(l) = zero
        ENDDO
        DO L=KD,K
          EVP(L) = zero
          ETZ(L) = zero
        ENDDO

      ENDIF

!
!===> CALCULATE GAMMAS  i.e. TENDENCIES PER UNIT CLOUD BASE MASSFLUX
!           Includes downdraft terms!

      avh = zero

!
!     Fraction of detrained condensate evaporated
!
!     tem1 = max(ZERO, min(HALF, (prl(kd)-FOUR_P2)*ONE_M2))
!     tem1 = max(ZERO, min(HALF, (prl(kd)-300.0)*0.005))
      tem1 = zero
!     tem1 = 1.0
!     if (kd1 == kbl) tem1 = 0.0
!
      tem2 = one - tem1
      TEM  = DET * QIL(KD)


      st1 = (HCC+ALHF*TEM-ETA(KD)*HST(KD)) / (one+gam(KD))
      DS  = ETA(KD1) * (HOS- HOL(KD)) - ALHL*(QOS - QOL(KD))
      DH  = ETA(KD1) * (HOS- HOL(KD))


      GMS(KD) = (DS + st1 - tem1*det*alhl-tem*alhf) * PRI(KD)
      GMH(KD) = PRI(KD) * (HCC-ETA(KD)*HOS + DH)
!
!      TENDENCY FOR SUSPENDED ENVIRONMENTAL ICE AND/OR LIQUID WATER
!
      QLL(KD) = (tem2*(DET-TEM) + ETA(KD1)*(QLS-CLL(KD))                &
     &        + (one-QIL(KD))*dlq(kd) - ETA(KD)*QLS ) * PRI(KD)

      QIL(KD) =     (tem2*TEM + ETA(KD1)*(QIS-CIL(KD))                  &
     &        + QIL(KD)*dlq(kd) - ETA(KD)*QIS ) * PRI(KD)
!
      GHD(KD) = zero
      GSD(KD) = zero
!
      DO L=KD1,K
         lm1 = l - 1
         ST1 = ONE - ALFINT(L,1)
         ST2 = ONE - ALFINT(L,2)
         ST3 = ONE - ALFINT(L,3)
         ST4 = ONE - ALFINT(L,4)
         ST5 = ONE - ALFIND(L)
         HB       = ALFINT(L,1)*HOL(LM1) + ST1*HOL(L)
         QB       = ALFINT(L,2)*QOL(LM1) + ST2*QOL(L)

         TEM      = ALFINT(L,4)*CIL(LM1) + ST4*CIL(L)
         TEM2     = ALFINT(L,3)*CLL(LM1) + ST3*CLL(L)
 
         TEM1     = ETA(L) * (TEM - CIL(L))
         TEM3     = ETA(L) * (TEM2 - CLL(L))

         HBD      = ALFIND(L)*HOL(LM1) + ST5*HOL(L)
         QBD      = ALFIND(L)*QOL(LM1) + ST5*QOL(L)

         TEM5     = ETD(L) * (HOD(L) - HBD)
         TEM6     = ETD(L) * (QOD(L) - QBD)
!
         DH       = ETA(L) * (HB - HOL(L)) + TEM5
         DS       = DH - ALHL * (ETA(L) * (QB - QOL(L)) + TEM6)

         GMH(L)   = DH * PRI(L)
         GMS(L)   = DS * PRI(L)
!
         GHD(L)   = TEM5 * PRI(L)
         GSD(L)   = (TEM5 - ALHL * TEM6) * PRI(L)
!
         QLL(L)   = (TEM3 + (one-QIL(L))*dlq(l)) * PRI(L)
         QIL(L)   = (TEM1 + QIL(L)*dlq(l)) * PRI(L)

         TEM1     = ETA(L) * (CIL(LM1) - TEM)
         TEM3     = ETA(L) * (CLL(LM1) - TEM2)

         DH       = ETA(L) * (HOL(LM1) - HB) - TEM5
         DS       = DH - ALHL * ETA(L) * (QOL(LM1) - QB)                &
     &                 + ALHL * (TEM6 - EVP(LM1))

         GMH(LM1) = GMH(LM1) + DH * PRI(LM1)
         GMS(LM1) = GMS(LM1) + DS * PRI(LM1)
!
         GHD(LM1) = GHD(LM1) - TEM5 * PRI(LM1)
         GSD(LM1) = GSD(LM1) - (TEM5-ALHL*(TEM6-EVP(LM1))) * PRI(LM1)

         QIL(LM1) = QIL(LM1) + TEM1 * PRI(LM1)
         QLL(LM1) = QLL(LM1) + TEM3 * PRI(LM1)
!
        avh = avh + gmh(lm1)*(prs(l)-prs(lm1))

      ENDDO
!
      HBD    = HOL(K)
      QBD    = QOL(K)
      TEM5   =  ETD(KP1) * (HOD(KP1) - HBD)
      TEM6   =  ETD(KP1) * (QOD(KP1) - QBD)
      DH     = - TEM5
      DS     = DH + ALHL * TEM6
      TEM1   = DH * PRI(K)
      TEM2   = (DS - ALHL * EVP(K)) * PRI(K)
      GMH(K) = GMH(K) + TEM1
      GMS(K) = GMS(K) + TEM2
      GHD(K) = GHD(K) + TEM1
      GSD(K) = GSD(K) + TEM2

!
      avh    = avh + gmh(K)*(prs(KP1)-prs(K))
!
      tem4   = - GRAVFAC * pris
      TX1    = DH * tem4
      TX2    = DS * tem4
!
      DO L=KBL,K
        GMH(L) = GMH(L) + TX1
        GMS(L) = GMS(L) + TX2
        GHD(L) = GHD(L) + TX1
        GSD(L) = GSD(L) + TX2
!
        avh = avh + tx1*(prs(l+1)-prs(l))
      ENDDO
!
!***********************************************************************
!***********************************************************************

!===>  KERNEL (AKM) CALCULATION BEGINS

!===>  MODIFY SOUNDING WITH UNIT MASS FLUX
!
      DO L=KD,K

         TEM1   = GMH(L)
         TEM2   = GMS(L)
         HOL(L) = HOL(L) +  TEM1*TESTMB
         QOL(L) = QOL(L) + (TEM1-TEM2)  * TESTMBOALHL
         HST(L) = HST(L) +  TEM2*(ONE+GAM(L))*TESTMB
         QST(L) = QST(L) +  TEM2*GAM(L) * TESTMBOALHL
         CLL(L) = CLL(L) + QLL(L) * TESTMB
         CIL(L) = CIL(L) + QIL(L) * TESTMB
      ENDDO
!
      if (alm > zero) then
        HOS     = HOS + GMH(KD)  * TESTMB
        QOS     = QOS + (GMH(KD)-GMS(KD)) * TESTMBOALHL
        QLS     = QLS + QLL(KD) * TESTMB
        QIS     = QIS + QIL(KD) * TESTMB
      else
        st2 = one - st1s
        HOS     = HOS + (st1s*GMH(KD)+st2*GMH(KD1))  * TESTMB
        QOS     = QOS + (st1s * (GMH(KD)-GMS(KD))                       &
     &                +  st2  * (GMH(KD1)-GMS(KD1))) * TESTMBOALHL
        HST(kd) = HST(kd) + (st1s*GMS(kd)*(ONE+GAM(kd))                 &
     &                    +  st2*gms(kd1)*(ONE+GAM(kd1))) * TESTMB
        QST(kd) = QST(kd) + (st1s*GMS(kd)*GAM(kd)                       &
     &                    +  st2*gms(kd1)*gam(kd1)) * TESTMBOALHL

        QLS     = QLS + (st1s*QLL(KD)+st2*QLL(KD1)) * TESTMB
        QIS     = QIS + (st1s*QIL(KD)+st2*QIL(KD1)) * TESTMB
      endif

!
      TEM = PRL(Kmaxp1) - PRL(Kmax)
      HBL = HOL(Kmax) * TEM
      QBL = QOL(Kmax) * TEM
      QLB = CLL(Kmax) * TEM
      QIB = CIL(Kmax) * TEM
      DO L=KmaxM1,KBL,-1
        TEM = PRL(L+1) - PRL(L)
        HBL = HBL + HOL(L) * TEM
        QBL = QBL + QOL(L) * TEM
        QLB = QLB + CLL(L) * TEM
        QIB = QIB + CIL(L) * TEM
      ENDDO
      HBL = HBL * PRISM
      QBL = QBL * PRISM
      QLB = QLB * PRISM
      QIB = QIB * PRISM

!     if (ctei .and. sgcs(kd) > 0.65) then
!        hbl = hbl * hpert_fac
!        qbl = qbl * hpert_fac
!     endif
 
!***********************************************************************

!===>  CLOUD WORKFUNCTION FOR MODIFIED SOUNDING, THEN KERNEL (AKM)
!
      AKM = ZERO
      TX1 = ZERO
      QTL = QST(KB1) - GAF(KB1)*HST(KB1)
      QTV = QBL
      HCC = HBL
      TX2 = HCC
      TX4 = (ALHF*half)*MAX(ZERO,MIN(ONE,(TCR-TCL-TOL(KB1))*TCRF))
!
      qtv = qbl
      tx1 = qib + qlb
!

      DO L=KB1,KD1,-1
         lm1 = l - 1
         lp1 = l + 1
         DEL_ETA = ETA(L) - ETA(LP1)
         HCCP = HCC + DEL_ETA*HOL(L)
!
         QTLP = QST(LM1) - GAF(LM1)*HST(LM1)
         QTVP = half * ((QTLP+QTL)*ETA(L) + (GAF(L)+GAF(LM1))*HCCP)

         DETP = (BKC(L)*TX1 - (QTVP-QTV)                                &
     &        +  DEL_ETA*(QOL(L)+CLL(L)+CIL(L))                         &
     &        +  ETA(L)*Q0U(L) + ETA(LP1)*Q0D(L)) * AKC(L)
         IF (DETP <= ZERO) cnvflg = .TRUE.

         ST1  = HST(L) - LTL(L)*NU*(QST(L)-QOL(L))

         TEM2 = (ALHF*half)*MAX(ZERO,MIN(ONE,(TCR-TCL-TOL(LM1))*TCRF))
         TEM1 = HCCP + DETP * (TEM2+TX4)

         ST2  = LTL(L) * VTF(L)
         TEM5 = CLL(L) + CIL(L)
         AKM  = AKM +                                                   &
     &     (  (TX2  -ETA(LP1)*ST1-ST2*(TX1-TEM5*eta(lp1))) * DLB(L)     &
     &      + (TEM1 -ETA(L  )*ST1-ST2*(DETP-TEM5*eta(l)))  * DLT(L) )
!
         HCC  = HCCP
         TX1  = DETP
         TX2  = TEM1
         QTL  = QTLP
         QTV  = QTVP
         TX4  = TEM2
      ENDDO
!
      if (cnvflg) return
!
!  Eventhough we ignore the change in lambda, we still assume
!  that the cLoud-top contribution is zero; as though we still
!  had non-bouyancy there.
!
!
      ST1  = HST(KD)  - LTL(KD)*NU*(QST(KD)-QOS)
      ST2  = LTL(KD)  * VTF(KD)
      TEM5 = (QLS + QIS) * eta(kd1)
      AKM  = AKM + HALF * (TX2-ETA(KD1)*ST1-ST2*(TX1-TEM5)) * DLB(KD)
!
      AKM  = (AKM - WFN) * TESTMBI


!***********************************************************************

!===>   MASS FLUX
!
      AMB  = - (WFN-ACR) / AKM
!
!===>   RELAXATION AND CLIPPING FACTORS
!
      AMB  = AMB * CLP * rel_fac

!!!   if (DDFT) AMB = MIN(AMB, ONE/CLDFRD)
       
!===>   SUB-CLOUD LAYER DEPTH LIMIT ON MASS FLUX

      AMBMAX = (PRL(KMAXP1)-PRL(KBL))*(FRACBL*GRAVCON)
      AMB    = MAX(MIN(AMB, AMBMAX),ZERO)

!***********************************************************************
!*************************RESULTS***************************************
!***********************************************************************

!===>  PRECIPITATION AND CLW DETRAINMENT
!
      if (amb > zero) then

!
!       if (wvl(kd) > zero) then
!         tx1 = one - amb * eta(kd) / (rho(kd)*wvl(kd))
!         sigf(kd) =  max(zero, min(one, tx1 * tx1))
!       endif
        if (aw_scal) then
          tx1 = (0.2_kp / max(alm, 1.0e-5_kp))
          tx2 = one - min(one, pi * tx1 * tx1 / area)

          tx2 = tx2 * tx2 

! comnet out the following for now - 07/23/18
!         do l=kd1,kbl
!           lp1 = min(K, l+1)
!           if (wvl(l) > zero .and. wvl(lp1) > zero) then
!             tx1     = one - amb *  (eta(l)+eta(lp1))
!    &                      / ((wvl(l)+wvl(lp1))*rho(l)*grav)
!             sigf(l) = max(zero, min(one, tx1 * tx1))
!           else
!             sigf(l) = min(one,tx2)
!           endif
!           sigf(l) = max(sigf(l), tx2)
!         enddo
!         sigf(kd) = sigf(kd1)
!         if (kbl < k) then
!           sigf(kbl+1:k) = sigf(kbl)
!         endif
          sigf(kd:k) = tx2
        else
          sigf(kd:k) = one
        endif

        tx1 = max(1.0e-6_kp, abs(gms(kd) * onebcp * sigf(kd)))
        amb = min(tx1*amb, tfrac_max*toi(kd)) / tx1

!
        avt = zero
        avq = zero
        avr = dof * sigf(kbl)
!
        DSFC = DSFC + AMB * ETD(K) * (one/DT) * sigf(kbl)
!
        DO L=K,KD,-1
          PCU(L) = PCU(L) + AMB*RNN(L)*sigf(l)      !  (A40)
          avr    = avr + rnn(l) * sigf(l)
        ENDDO
        pcu(k) = pcu(k) + amb * dof * sigf(kbl)
!
!===> TEMPARATURE AND Q CHANGE AND CLOUD MASS FLUX DUE TO CLOUD TYPE KD
!
        TX1 = AMB * ONEBCP
        TX2 = AMB * ONEOALHL
        DO L=KD,K
          delp   = prs(l+1) - prs(l)
          tx3    = amb * sigf(l)
          ST1    = GMS(L) * TX1 * sigf(l)
          TOI(L) = TOI(L) + ST1
          TCU(L) = TCU(L) + ST1
          TCD(L) = TCD(L) + GSD(L) * TX1 * sigf(l)
!
          st1 = st1 - ELOCP * (QIL(L) + QLL(L)) * tx3

          avt = avt + st1 * delp

          FLX(L)  = FLX(L)  + ETA(L) * tx3
          FLXD(L) = FLXD(L) + ETD(L) * tx3
!
          QII(L)  = QII(L) + QIL(L) * tx3
          TEM     = zero

          QLI(L)  = QLI(L) + QLL(L) * tx3 + TEM

          ST1     = (GMH(L)-GMS(L)) * TX2 * sigf(l)

          QOI(L)  = QOI(L) + ST1
          QCU(L)  = QCU(L) + ST1
          QCD(L)  = QCD(L) + (GHD(L)-GSD(L)) * TX2 * sigf(l)
!
          avq = avq + (st1 + (QLL(L)+QIL(L))*tx3) * delp
!         avq = avq + st1 * (prs(l+1)-prs(l))
!         avr = avr + (QLL(L) + QIL(L)*(1+alhf/alhl))
          avr = avr + (QLL(L) + QIL(L)) * delp * sigf(l) * gravcon

!      Correction for negative condensate!
          if (qii(l) < zero) then
            tem    = qii(l) * elfocp
            QOI(L) = QOI(L) + qii(l)
            qcu(l) = qcu(l) + qii(l)
            toi(l) = toi(l) - tem
            tcu(l) = tcu(l) - tem
            qii(l) = zero
          endif
          if (qli(l) < zero) then
            tem    = qli(l) * elocp
            QOI(L) = QOI(L) + qli(l)
            qcu(l) = qcu(l) + qli(l)
            toi(l) = toi(l) - tem
            tcu(l) = tcu(l) - tem
            qli(l) = zero
          endif

        ENDDO
        avr = avr * amb
!
!      Correction for negative condensate!
!     if (advcld) then
!       do l=kd,k
!         if (qli(l) < zero) then
!           qoi(l) = qoi(l) + qli(l)
!           toi(l) = toi(l) - (alhl/cp) * qli(l)
!           qli(l) = zero
!         endif
!         if (qii(l) < zero) then
!           qoi(l) = qoi(l) + qii(l)
!           toi(l) = toi(l) - ((alhl+alhf)/cp) * qii(l)
!           qii(l) = zero
!         endif
!       enddo
!     endif

!
!
        TX1 = zero
        TX2 = zero
!
        IF (REVAP) THEN !     REEVAPORATION OF FALLING CONVECTIVE RAIN
!
         tem = zero
         do l=kd,kbl
           IF (L < IDH .or. (.not. DDFT)) THEN
             tem = tem + amb * rnn(l) * sigf(l)
           endif
         enddo
         tem = tem + amb * dof * sigf(kbl)
         tem = tem * (3600.0_kp/dt)
         tem1 = sqrt(max(one, min(100.0_kp,(6.25e10_kp/max(area,one))))) ! 20110530

         clfrac = max(ZERO, min(half, rknob*clf(tem)*tem1))
         cldfrd = clfrac
!
         DO L=KD,KBL         ! Testing on 20070926
!                                                 for L=KD,K
           IF (L >= IDH .AND. DDFT) THEN
             tem    = amb * sigf(l)
             TX2    = TX2 + tem * RNN(L)
             CLDFRD = MIN(tem*CLDFR(L), clfrac)
           ELSE
             TX1 = TX1 + AMB * RNN(L) * sigf(l)
           ENDIF
           tx4 = zfac * phil(l)
           tx4 = (one - tx4 * (one - half*tx4)) * afc
!
           IF (TX1 > zero .OR. TX2 > zero) THEN
             TEQ     = TOI(L)
             QEQ     = QOI(L)
             PL      = half * (PRL(L+1)+PRL(L))

             ST1     = MAX(ZERO, MIN(ONE, (TCR-TEQ)*TCRF))
             ST2     = ST1*ELFOCP + (one-ST1)*ELOCP

             CALL QSATCN ( TEQ,PL,QSTEQ,DQDT)
!
             DELTAQ = half * (QSTEQ*rhc_ls(l)-QEQ) / (one+ST2*DQDT)
!
             QEQ    = QEQ + DELTAQ
             TEQ    = TEQ - DELTAQ*ST2
!
             TEM1   = MAX(ZERO, MIN(ONE, (TCR-TEQ)*TCRF))
             TEM2   = TEM1*ELFOCP + (one-TEM1)*ELOCP

             CALL QSATCN ( TEQ,PL,QSTEQ,DQDT)
!
             DELTAQ = (QSTEQ*rhc_ls(l)-QEQ) / (one+TEM2*DQDT)
!
             QEQ    = QEQ + DELTAQ
             TEQ    = TEQ - DELTAQ*TEM2

             IF (QEQ > QOI(L)) THEN
               POTEVAP = (QEQ-QOI(L))*(PRL(L+1)-PRL(L))*GRAVCON

               tem4    = zero
               if (tx1 > zero)                                          &
     &         TEM4    = POTEVAP * (one - EXP( tx4*TX1**0.57777778_kp ))
               ACTEVAP = MIN(TX1, TEM4*CLFRAC)

               if (tx1 < rainmin*dt) actevap = min(tx1, potevap)
!
               tem4    = zero
               if (tx2 > zero)                                          &
     &         TEM4    = POTEVAP * (one - EXP( tx4*TX2**0.57777778_kp ))
               TEM4    = min(MIN(TX2, TEM4*CLDFRD), potevap-actevap)
               if (tx2 < rainmin*dt) tem4 = min(tx2, potevap-actevap)
!
               TX1     = TX1 - ACTEVAP
               TX2     = TX2 - TEM4
               ST1     = (ACTEVAP+TEM4) * PRI(L)
               QOI(L)  = QOI(L) + ST1
               QCU(L)  = QCU(L) + ST1
!

               ST1     = ST1 * ELOCP
               TOI(L)  = TOI(L) - ST1 
               TCU(L)  = TCU(L) - ST1
             ENDIF
           ENDIF
         ENDDO
!
          CUP = CUP + TX1 + TX2 + DOF * AMB * sigf(kbl)
        ELSE
          DO L=KD,K
            TX1 = TX1 + AMB * RNN(L) * sigf(l)
          ENDDO
          CUP = CUP + TX1 + DOF * AMB * sigf(kbl)
        ENDIF
!
!    Convective transport (mixing) of passive tracers
!
        if (NTRC > 0) then
          do l=kd,km1
            if (etz(l) /= zero) etzi(l) = one / etz(l)
          enddo
          DO N=1,NTRC        ! Tracer loop ; first two are u and v

            DO L=KD,K
              HOL(L) = ROI(L,N)
            ENDDO
!
            HCC     = RBL(N)
            HOD(KD) = HOL(KD)
!      Compute downdraft properties for the tracer
            DO L=KD1,K
              lm1 = l - 1
              ST1 = ONE - ALFIND(L)
              HB  = ALFIND(L)  * HOL(LM1) + ST1 * HOL(L)
              IF (ETZ(LM1) /= ZERO) THEN
                TEM = ETZI(LM1)
                IF (ETD(L)  > ETD(LM1)) THEN
                 HOD(L) = (ETD(LM1)*(HOD(LM1)-HOL(LM1))                 &
     &                   +  ETD(L)  *(HOL(LM1)-HB) +  ETZ(LM1)*HB) * TEM
                ELSE
                 HOD(L) = (ETD(LM1)*(HOD(LM1)-HB) + ETZ(LM1)*HB) * TEM
                ENDIF
               ELSE
                HOD(L) = HB
              ENDIF
            ENDDO
             
            DO L=KB1,KD,-1
              HCC = HCC + (ETA(L)-ETA(L+1))*HOL(L)
            ENDDO
!
!         Scavenging -- fscav   - fraction scavenged [km-1]
!                       delz    - distance from the entrainment to detrainment layer [km]
!                       fnoscav - the fraction not scavenged
!                                 following Liu et al. [JGR,2001] Eq 1

            if (FSCAV_(N) > zero) then
              DELZKM = ( PHIL(KD) - PHIH(KD1) ) *(onebg*0.001_kp)
              FNOSCAV = exp(- FSCAV_(N) * DELZKM)
            else
              FNOSCAV = one
            endif

            GMH(KD) = PRI(KD) * (HCC-ETA(KD)*HOL(KD)) * trcfac(kd,n)    &
     &                                                * FNOSCAV
            DO L=KD1,K
              if (FSCAV_(N) > zero) then
                DELZKM = ( PHIL(KD) - PHIH(L+1) ) *(onebg*0.001_kp)
                FNOSCAV = exp(- FSCAV_(N) * DELZKM)
              endif
              lm1      = l - 1
              ST1      = ONE - ALFINT(L,N+4)
              ST2      = ONE - ALFIND(L)
              HB       = ALFINT(L,N+4) * HOL(LM1) + ST1 * HOL(L)
              HBD      = ALFIND(L) * HOL(LM1) + ST2 * HOL(L)
              TEM5     = ETD(L)    * (HOD(L) - HBD)
              DH       = ETA(L)    * (HB - HOL(L))   * FNOSCAV + TEM5
              GMH(L  ) = DH * PRI(L) * trcfac(l,n)
              DH       = ETA(L)    * (HOL(LM1) - HB) * FNOSCAV - TEM5
              GMH(LM1) = GMH(LM1)  + DH * PRI(LM1) * trcfac(l,n)
            ENDDO
!
            st2 = zero
            DO L=KD,K
              ST1 = GMH(L)*AMB*sigf(l) + st2
              st3 = HOL(L) + st1
              st2 = st3 - trcmin(n) ! if trcmin is defined limit change
              if (st2 < zero) then
                ROI(L,N) = trcmin(n)
                RCU(L,N) = RCU(L,N) + ST1
                if (l < k)                                              &
     &          st2 = st2 * (prl(l+1)-prl(l))*pri(l+1) * (cmb2pa/grav)
              else
                ROI(L,N) = ST3
                RCU(L,N) = RCU(L,N) + ST1
                st2 = zero
              endif
              
            ENDDO
          ENDDO                             ! Tracer loop NTRC
        endif
      endif             ! amb > zero

      RETURN
      end subroutine cloud

!>\ingroup rascnv_schm
      SUBROUTINE DDRFT(                                                 &
     &                  K,   KP1, KD                                    &
     &,                 TLA, ALFIND, wcbase                             &
     &,                 TOL, QOL, HOL, PRL, QST, HST, GAM, GAF          &
!    &,                 TOL, QOL, HOL, PRL, QST, HST, GAM, GAF, HBL, QBL&
     &,                 QRB, QRT, BUY, KBL, IDH, ETA, RNN, ETAI         &
     &,                 ALM, WFN, TRAIN, DDFT                           &
     &,                 ETD, HOD, QOD, EVP, DOF, CLDFRD, WCB            &
     &,                 GMS, GSD, GHD, wvlu)                   

!
!***********************************************************************
!******************** Cumulus Downdraft Subroutine *********************
!****************** Based on Cheng and Arakawa (1997)  ****** **********
!************************ SUBROUTINE DDRFT  ****************************
!*************************  October 2004  ******************************
!***********************************************************************
!***********************************************************************
!************* Shrinivas.Moorthi@noaa.gov (301) 683-3718 ***************
!***********************************************************************
!***********************************************************************
!23456789012345678901234567890123456789012345678901234567890123456789012
!
!===>  TOL(K)     INPUT   TEMPERATURE            KELVIN
!===>  QOL(K)     INPUT   SPECIFIC HUMIDITY      NON-DIMENSIONAL

!===>  PRL(KP1)   INPUT   PRESSURE @ EDGES       MB

!===>  K     INPUT   THE RISE & THE INDEX OF THE SUBCLOUD LAYER
!===>  KD    INPUT   DETRAINMENT LEVEL ( 1<= KD < K )          
!     
      IMPLICIT NONE
!
!  INPUT ARGUMENTS
!
      INTEGER K, KP1, KD, KBL
      real(kind=kind_phys) ALFIND(K), wcbase

      real(kind=kind_phys), dimension(kd:k)   ::  HOL, QOL, HST, QST    &
     &,                                           TOL, QRB, QRT, RNN    &
     &,                                           RNS, ETAI
      real(kind=kind_phys), dimension(kd:kp1) ::  GAF, BUY, GAM, ETA    &
     &,                                           PRL
!
!     real(kind=kind_phys)    HBL,     QBL,        PRIS                 &
!    &,                       TRAIN,   WFN,        ALM
!
!     TEMPORARY WORK SPACE
!
      real(kind=kind_phys), dimension(KD:K)   :: RNF, WCB,  EVP, STLT   &
     &,                                          GHD, GSD,  CLDFRD      &
     &,                                          GQW, QRPI, QRPS, BUD

      real(kind=kind_phys), dimension(KD:KP1) :: QRP, WVL, WVLU, ETD    &
     &,                                          HOD, QOD, ROR,  GMS

      real(kind=kind_phys) TL,     PL,     QL,      QS,   DQS,  ST1     &
     &,                    QQQ,            DEL_ETA, HB,   QB,   TB      &
     &,                    TEM,    TEM1,   TEM2,    TEM3, TEM4, ST2     &
     &,                    ERRMIN, ERRMI2, ERRH,    ERRW, ERRE, TEM5    &
     &,                    TEM6,   HBD,    QBD,     TX1,  TX2,  TX3     &
     &,                    TX4,    TX5,    TX6,     TX7,  TX8,  TX9     &
     &,                    WFN,    ALM,             AL2                 &
     &,                    TRAIN,  GMF,    ONPG,    CTLA, VTRM          &
     &,                    RPART,  QRMIN,  AA1,     BB1,  CC1,   DD1    &
!    &,                    WC2MIN, WCMIN,  WCBASE,  F2,   F3,    F5     &
     &,                    WC2MIN, WCMIN,           F2,   F3,    F5     &
     &,                    GMF1,   GMF5,   QRAF,    QRBF, del_tla       &
     &,                    TLA,    STLA,   CTL2,    CTL3                &
!    &,                    TLA,    STLA,   CTL2,    CTL3, ASIN          &
!    &,                    RNT,    RNB,    ERRQ,    RNTP, QRPF,  VTPF   &
     &,                    RNT,    RNB,    ERRQ,    RNTP                &
     &,                    EDZ,    DDZ,    CE,      QHS,  FAC,   FACG   &
     &,                    RSUM1,  RSUM2,  RSUM3,   CEE,  DOF,   DOFW
!    &,                    sialf

      INTEGER              I, L,  N, IX, KD1, II, kb1, IP1, JJ, ntla    &
     &,                    IT, KM1, KTEM, KK, KK1, LM1, LL, LP1         &
     &,                    IDW, IDH, IDN(K), idnm, itr
!
      parameter (ERRMIN=0.0001_kp, ERRMI2=0.1_kp*ERRMIN)
!     parameter (ERRMIN=0.00001, ERRMI2=0.1*ERRMIN)
!
!     real (kind=kind_phys), parameter :: PIINV=one/PI, pio2=half*pi
!
      parameter (ONPG=one+half, GMF=one/ONPG, RPART=zero)
!     parameter (ONPG=1.0+0.5, GMF=1.0/ONPG, RPART=1.0)
!     parameter (ONPG=1.0+0.5, GMF=1.0/ONPG, RPART=0.5)
!     PARAMETER (AA1=1.0, BB1=1.5, CC1=1.1, DD1=0.85, F3=CC1, F5=2.5)
!     PARAMETER (AA1=2.0, BB1=1.5, CC1=1.1, DD1=0.85, F3=CC1, F5=2.5)
      PARAMETER (AA1=1.0_kp, BB1=1.0_kp, CC1=1.0_kp, DD1=1.0_kp,        &
     &           F3=CC1,  F5=1.0_kp)
      parameter (QRMIN=1.0e-6_kp, WC2MIN=0.01_kp, GMF1=GMF/AA1, GMF5=GMF/F5)
!     parameter (QRMIN=1.0E-6, WC2MIN=1.00, GMF1=GMF/AA1, GMF5=GMF/F5)
      parameter (WCMIN=sqrt(wc2min))
!     parameter (sialf=0.5)
!
      integer, parameter :: itrmu=25,  itrmd=25                         &
     &,                     itrmin=15, itrmnd=12, numtla=2

!     uncentering for vvel in dd
      real(kind=kind_phys), parameter :: ddunc1=0.25_kp                 &
     &,                                  ddunc2=one-ddunc1              &
!    &,                                  ddunc1=0.4,  ddunc2=one-ddunc1 &
!    &,                                  ddunc1=0.3,  ddunc2=one-ddunc1 &
     &,                      VTPEXP=-0.3636_kp                          &
     &,                      VTP=36.34_kp*SQRT(1.2_kp)*(0.001_kp)**0.1364_kp
!
!     real(kind=kind_phys) EM(K*K), ELM(K)
      real(kind=kind_phys) ELM(K), AA(KD:K,KD:KP1), QW(KD:K,KD:K)       &
     &,                    VT(2),  VRW(2), TRW(2), QA(3), WA(3)

      LOGICAL SKPUP, cnvflg, DDFT, UPDRET, DDLGK

!***********************************************************************


      KD1    = KD + 1
      KM1    = K  - 1
      KB1    = KBL - 1
!
!     VTP    = 36.34*SQRT(1.2)* (0.001)**0.1364
!     VTPEXP = -0.3636
!     PIINV  = 1.0 / PI
!     PICON  = PIO2 * ONEBG
!
!     Compute Rain Water Budget of the Updraft (Cheng and Arakawa, 1997)
!
      CLDFRD = zero
      RNTP   = zero
      DOF    = zero
      ERRQ   = 10.0_kp
      RNB    = zero
      RNT    = zero
      TX2    = PRL(KBL)
!
      TX1      = (PRL(KD) + PRL(KD1)) * half
      ROR(KD)  = CMPOR*TX1 / (TOL(KD)*(one+NU*QOL(KD)))
!     GMS(KD)  = VTP * ROR(KD) ** VTPEXP
      GMS(KD)  = VTP * VTPF(ROR(KD))
!
      QRP(KD)  = QRMIN
!
      TEM      = TOL(K) * (one + NU * QOL(K))
      ROR(KP1) = half * CMPOR * (PRL(KP1)+PRL(K)) / TEM
      GMS(KP1) = VTP * VTPF(ROR(KP1))
      QRP(KP1) = QRMIN
!
      kk = kbl
      DO L=KD1,K
        TEM = half * (TOL(L)+TOL(L-1))                                  &
     &      * (one + (half*NU) * (QOL(L)+QOL(L-1)))
        ROR(L) = CMPOR * PRL(L) / TEM
!       GMS(L) = VTP * ROR(L) ** VTPEXP
        GMS(L) = VTP * VTPF(ROR(L))
        QRP(L) = QRMIN
        if (buy(l) <= zero .and. kk == KBL) then
          kk = l
        endif
      ENDDO
      if (kk /= kbl) then
        do l=kk,kbl
          buy(l) = 0.9_kp * buy(l-1)
        enddo
      endif
!
      do l=kd,k
        qrpi(l) = buy(l)
      enddo
      do l=kd1,kb1
        buy(l) = 0.25_kp * (qrpi(l-1)+qrpi(l)+qrpi(l)+qrpi(l+1))
      enddo
      
!
!     CALL ANGRAD(TX1, ALM, STLA, CTL2, AL2, PI, TLA, TX2, WFN, TX3)
      tx1 = 1000.0_kp + tx1 - prl(kp1)
!     CALL ANGRAD(TX1, ALM,  AL2, TLA, TX2, WFN, TX3)
      CALL ANGRAD(TX1, ALM,  AL2, TLA)
!
!    Following Ucla approach for rain profile
!
      F2      = (BB1+BB1)*ONEBG/(PI*0.2_kp)
!     WCMIN   = SQRT(WC2MIN)
!     WCBASE  = WCMIN
!
!     del_tla = TLA * 0.2
!     del_tla = TLA * 0.25
      del_tla = TLA * 0.3_kp
      TLA     = TLA - DEL_TLA
!
      DO L=KD,K
        RNF(L)   = zero
        RNS(L)   = zero
        STLT(L)  = zero
        GQW(L)   = zero
        QRP(L)   = QRMIN
        DO N=KD,K
          QW(N,L) = zero
        ENDDO
      ENDDO
!     DO L=KD,KP1
!       WVL(L)   = zero
!     ENDDO
!
!-----QW(N,L) = D(W(N)*W(N))/DQR(L)
!
      KK = KBL
      QW(KD,KD) = -QRB(KD)  * GMF1
      GHD(KD)   = ETA(KD)   * ETA(KD)
      GQW(KD)   = QW(KD,KD) * GHD(KD)
      GSD(KD)   = ETAI(KD)  * ETAI(KD)
!
      GQW(KK)   = - QRB(KK-1) * (GMF1+GMF1)
!
      WCB(KK)   = WCBASE * WCBASE

      TX1       = WCB(KK)
      GSD(KK)   = one
      GHD(KK)   = one
!
      TEM       = GMF1 + GMF1
      DO L=KB1,KD1,-1
        GHD(L)  = ETA(L)  * ETA(L)
        GSD(L)  = ETAI(L) * ETAI(L)
        GQW(L)  = - GHD(L) * (QRB(L-1)+QRT(L)) * TEM
        QW(L,L) = - QRT(L) * TEM
!
        st1     = half * (eta(l) + eta(l+1))
        TX1     = TX1 + BUY(L) * TEM * (qrb(l)+qrt(l)) * st1 * st1
        WCB(L)  = TX1 * GSD(L)
      ENDDO
!
      TEM1        = (QRB(KD) + QRT(KD1) + QRT(KD1)) * GMF1
      GQW(KD1)    = - GHD(KD1) * TEM1
      QW(KD1,KD1) = - QRT(KD1) * TEM
      st1         = half * (eta(kd) + eta(kd1))
      WCB(KD)     = (TX1 + BUY(KD)*TEM*qrb(kd)*st1*st1) * GSD(KD)
!
      DO L=KD1,KBL
        DO N=KD,L-1
          QW(N,L) = GQW(L) * GSD(N)
        ENDDO
      ENDDO
      QW(KBL,KBL) = zero
!
      do ntla=1,numtla    ! numtla is the the maximimu number of tilting angle tries
                          ! ------
!       if (errq < 1.0 .or. tla > 45.0) cycle
        if (errq < 0.1_kp .or. tla > 45.0_kp) cycle
!
        tla  = tla + del_tla
        STLA = SIN(TLA*deg2rad)     ! sine of tilting angle
        CTL2 = one - STLA * STLA    ! cosine square of tilting angle
!
        STLA = F2        * STLA * AL2
        CTL2 = DD1       * CTL2
        CTL3 = 0.1364_kp * CTL2
!
        DO L=KD,K
          RNF(L)  = zero
          STLT(L) = zero
          QRP(L)  = QRMIN
        ENDDO
        DO L=KD,KP1
          WVL(L)  = zero
        ENDDO
        WVL(KBL)  = WCBASE
        STLT(KBL) = one / WCBASE
!
        DO L=KD,KP1
          DO N=KD,K
            AA(N,L) = zero
          ENDDO
        ENDDO
!
        SKPUP = .FALSE.
!
        DO ITR=1,ITRMU               ! Rain Profile Iteration starts!
          IF (.NOT. SKPUP) THEN
!            wvlu = wvl
!
!-----CALCULATING THE VERTICAL VELOCITY
!
            TX1      = zero
            QRPI(KBL) = one / QRP(KBL)
            DO L=KB1,KD,-1
              TX1 = TX1    + QRP(L+1)*GQW(L+1)
              ST1 = WCB(L) + QW(L,L)*QRP(L) + TX1*GSD(L)
!             if (st1 > wc2min) then
              if (st1 > zero) then
                WVL(L) = max(ddunc1*SQRT(ST1) + ddunc2*WVL(L), wcmin)
!               WVL(L) = SQRT(ST1)
!               WVL(L) = max(half * (SQRT(ST1) + WVL(L)), wcmin)
!               qrp(l) = half*((wvl(l)*wvl(l)-wcb(l)-tx1*gsd(l))/qw(l,l)&
!    &                      + qrp(l))
              else

!               wvl(l) = 0.5*(wcmin+wvl(l))
!               wvl(l) = max(half*(wvl(l) + wvl(l+1)), wcmin)
                wvl(l) = max(wvl(l),wcmin)
                qrp(l) = (wvl(l)*wvl(l) - wcb(l) - tx1*gsd(l))/qw(l,l)
!               qrp(l) = half*((wvl(l)*wvl(l)-wcb(l)-tx1*gsd(l))/qw(l,l)&
!    &                      + qrp(l))
              endif
              qrp(l) = max(qrp(l), qrmin)

              STLT(L) = one / WVL(L)
              QRPI(L) = one / QRP(L)
            ENDDO
!
!-----CALCULATING TRW, VRW AND OF
!
!           VT(1)   = GMS(KD) * QRP(KD)**0.1364
            VT(1)   = GMS(KD) * QRPF(QRP(KD))
            TRW(1)  = ETA(KD) * QRP(KD) * STLT(KD)
            TX6     = TRW(1) * VT(1)
            VRW(1)  = F3*WVL(KD) - CTL2*VT(1)
            BUD(KD) = STLA * TX6 * QRB(KD) * half
            RNF(KD) = BUD(KD)
            DOF     = 1.1364_kp * BUD(KD) * QRPI(KD)
            DOFW    = -BUD(KD)  * STLT(KD)
!
            RNT     = TRW(1) * VRW(1)
            TX2     = zero
            TX4     = zero
            RNB     = RNT
            TX1     = half
            TX8     = zero
!
            IF (RNT >= zero) THEN
              TX3 = (RNT-CTL3*TX6) * QRPI(KD)
              TX5 = CTL2 * TX6 * STLT(KD)
            ELSE
              TX3 = zero
              TX5 = zero
              RNT = zero
              RNB = zero
            ENDIF
!
            DO L=KD1,KB1
              KTEM    = MAX(L-2, KD)
              LL      = L - 1
! 
!             VT(2)   = GMS(L) * QRP(L)**0.1364
              VT(2)   = GMS(L) * QRPF(QRP(L))
              TRW(2)  = ETA(L) * QRP(L) * STLT(L)
              VRW(2)  = F3*WVL(L) - CTL2*VT(2)
              QQQ     = STLA * TRW(2) * VT(2)
              ST1     = TX1  * QRB(LL)
              BUD(L)  = QQQ * (ST1 + QRT(L))
!
              QA(2)   = DOF
              WA(2)   = DOFW
              DOF     = 1.1364_kp * BUD(L) * QRPI(L)
              DOFW    = -BUD(L)   * STLT(L)
!
              RNF(LL) = RNF(LL) + QQQ * ST1
              RNF(L)  =           QQQ * QRT(L)
!
              TEM3    = VRW(1) + VRW(2)
              TEM4    = TRW(1) + TRW(2)
!
              TX6     = pt25 * TEM3 * TEM4
              TEM4    = TEM4 * CTL3
!
!-----BY QR ABOVE
!
!             TEM1    = pt25*(TRW(1)*TEM3 - TEM4*VT(1))*TX7
              TEM1    = pt25*(TRW(1)*TEM3 - TEM4*VT(1))*QRPI(LL)
              ST1     = pt25*(TRW(1)*(CTL2*VT(1)-VRW(2))                &
     &                      * STLT(LL) + F3*TRW(2))
!-----BY QR BELOW
              TEM2    = pt25*(TRW(2)*TEM3 - TEM4*VT(2))*QRPI(L)
              ST2     = pt25*(TRW(2)*(CTL2*VT(2)-VRW(1))                &
     &                      * STLT(L)  + F3*TRW(1))
!
!      From top to  the KBL-2 layer
!
              QA(1)   = TX2
              QA(2)   = QA(2) + TX3 - TEM1
              QA(3)   = -TEM2
!
              WA(1)   = TX4
              WA(2)   = WA(2) + TX5 - ST1
              WA(3)   = -ST2
!
              TX2     = TEM1
              TX3     = TEM2
              TX4     = ST1
              TX5     = ST2
!
              VT(1)   = VT(2)
              TRW(1)  = TRW(2)
              VRW(1)  = VRW(2)
!
              IF (WVL(KTEM) == WCMIN) WA(1) = zero
              IF (WVL(LL)   == WCMIN) WA(2) = zero
              IF (WVL(L)    == WCMIN) WA(3) = zero
              DO N=KTEM,KBL
                AA(LL,N) = (WA(1)*QW(KTEM,N) * STLT(KTEM)               &
     &                   +  WA(2)*QW(LL,N)   * STLT(LL)                 &
     &                   +  WA(3)*QW(L,N)    * STLT(L) ) * half
              ENDDO
              AA(LL,KTEM) = AA(LL,KTEM) + QA(1)
              AA(LL,LL)   = AA(LL,LL)   + QA(2)
              AA(LL,L)    = AA(LL,L)    + QA(3)
              BUD(LL)     = (TX8 + RNN(LL)) * half                      &
     &                      - RNB + TX6 - BUD(LL)
              AA(LL,KBL+1) = BUD(LL)
              RNB = TX6
              TX1 = one
              TX8 = RNN(LL)
            ENDDO
            L  = KBL
            LL = L - 1
!           VT(2)   = GMS(L) * QRP(L)**0.1364
            VT(2)   = GMS(L) * QRPF(QRP(L))
            TRW(2)  = ETA(L) * QRP(L) * STLT(L)
            VRW(2)  = F3*WVL(L) - CTL2*VT(2)
            ST1     = STLA * TRW(2) * VT(2) * QRB(LL)
            BUD(L)  = ST1

            QA(2)   = DOF
            WA(2)   = DOFW
            DOF     = 1.1364_kp * BUD(L) * QRPI(L)
            DOFW    = -BUD(L)   * STLT(L)
!
            RNF(LL) = RNF(LL) + ST1
!
            TEM3    = VRW(1) + VRW(2)
            TEM4    = TRW(1) + TRW(2)
!
            TX6     = pt25 * TEM3 * TEM4
            TEM4    = TEM4 * CTL3
!
!-----BY QR ABOVE
!
            TEM1    = pt25*(TRW(1)*TEM3 - TEM4*VT(1))*QRPI(LL)
            ST1     = pt25*(TRW(1)*(CTL2*VT(1)-VRW(2))                  &
     &                    * STLT(LL) + F3*TRW(2))
!-----BY QR BELOW
            TEM2    = pt25*(TRW(2)*TEM3 - TEM4*VT(2))*QRPI(L)
            ST2     = pt25*(TRW(2)*(CTL2*VT(2)-VRW(1))                  &
     &                    * STLT(L)  + F3*TRW(1))
!
!      For the layer next to the top of the boundary layer
!
            QA(1)   = TX2
            QA(2)   = QA(2) + TX3 - TEM1
            QA(3)   = -TEM2
!
            WA(1)   = TX4
            WA(2)   = WA(2) + TX5 - ST1
            WA(3)   = -ST2
!
            TX2     = TEM1
            TX3     = TEM2
            TX4     = ST1
            TX5     = ST2
!
            IDW     = MAX(L-2, KD)
!
            IF (WVL(IDW) == WCMIN) WA(1) = zero
            IF (WVL(LL)  == WCMIN) WA(2) = zero
            IF (WVL(L)   == WCMIN) WA(3) = zero
!
            KK = IDW
            DO N=KK,L
              AA(LL,N) = (WA(1)*QW(KK,N) * STLT(KK)                     &
     &                 +  WA(2)*QW(LL,N) * STLT(LL)                     &
     &                 +  WA(3)*QW(L,N)  * STLT(L) ) * half

            ENDDO
!
            AA(LL,IDW) = AA(LL,IDW) + QA(1)
            AA(LL,LL)  = AA(LL,LL)  + QA(2)
            AA(LL,L)   = AA(LL,L)   + QA(3)
            BUD(LL)    = (TX8+RNN(LL)) * half - RNB + TX6 - BUD(LL)
!
            AA(LL,L+1) = BUD(LL)
!
            RNB        = TRW(2) * VRW(2)
!
!      For the top of the boundary layer
!
            IF (RNB < zero) THEN
               KK    = KBL
               TEM   = VT(2) * TRW(2)
               QA(2) = (RNB - CTL3*TEM) * QRPI(KK)
               WA(2) = CTL2 * TEM * STLT(KK)
            ELSE
               RNB   = zero
               QA(2) = zero
               WA(2) = zero
            ENDIF
!
            QA(1) = TX2
            QA(2) = DOF + TX3 - QA(2)
            QA(3) = zero
!
            WA(1) = TX4
            WA(2) = DOFW + TX5 - WA(2)
            WA(3) = zero
!
            KK = KBL
            IF (WVL(KK-1) == WCMIN) WA(1) = zero
            IF (WVL(KK)   == WCMIN) WA(2) = zero
!
            DO II=1,2
               N = KK + II - 2
               AA(KK,N) = (WA(1)*QW(KK-1,N) * STLT(KK-1)                &
     &                  +  WA(2)*QW(KK,N)   * STLT(KK)) * half
            ENDDO
            FAC = half
            LL  = KBL
            L   = LL + 1
            LM1 = LL - 1
            AA(LL,LM1)  = AA(LL,LM1) + QA(1)
            AA(LL,LL)   = AA(LL,LL)  + QA(2)
            BUD(LL)     = half*RNN(LM1) - TX6 + RNB - BUD(LL)
            AA(LL,LL+1) = BUD(LL)
!
!-----SOLVING THE BUDGET EQUATIONS FOR DQR
!
            DO L=KD1,KBL
              LM1  = L - 1
              cnvflg = ABS(AA(LM1,LM1)) < ABS(AA(L,LM1))
              DO  N=LM1,KBL+1
                 IF (cnvflg) THEN
                    TX1       = AA(LM1,N)
                    AA(LM1,N) = AA(L,N)
                    AA(L,N)   = TX1
                 ENDIF
              ENDDO
              TX1 = AA(L,LM1) / AA(LM1,LM1)
              DO  N=L,KBL+1
                AA(L,N) = AA(L,N) - TX1 * AA(LM1,N)
              ENDDO
            ENDDO     
!
!-----BACK SUBSTITUTION AND CHECK IF THE SOLUTION CONVERGES
!
            KK = KBL
            KK1 = KK + 1
            AA(KK,KK1) = AA(KK,KK1) / AA(KK,KK)      !   Qr correction !
            TX2        = ABS(AA(KK,KK1)) * QRPI(KK)  !   Error Measure !
!
            KK = KBL + 1
            DO L=KB1,KD,-1
               LP1   = L + 1
               TX1  = zero
               DO N=LP1,KBL
                 TX1  = TX1 + AA(L,N) * AA(N,KK)
               ENDDO
               AA(L,KK) = (AA(L,KK) - TX1) / AA(L,L)       ! Qr correction !
               TX2      = MAX(TX2, ABS(AA(L,KK))*QRPI(L))  ! Error Measure !
            ENDDO
!
!           tem = 0.5
            if (tx2 > one .and. abs(errq-tx2) > 0.1_kp) then
              tem = half
!!          elseif (tx2 < 0.1) then
!!            tem = 1.2
            else
              tem = one
            endif
!
            DO L=KD,KBL
!              QRP(L) = MAX(QRP(L)+AA(L,KBL+1), QRMIN)
               QRP(L) = MAX(QRP(L)+AA(L,KBL+1)*tem, QRMIN)
            ENDDO
!
            IF (ITR < ITRMIN) THEN
               TEM = ABS(ERRQ-TX2) 
               IF (TEM >= ERRMI2 .AND. TX2 >= ERRMIN) THEN 
                 ERRQ  = TX2                              ! Further iteration !
               ELSE 
                 SKPUP = .TRUE.                           ! Converges      !
                 ERRQ  = zero                             ! Rain profile exists!
               ENDIF 
            ELSE
               TEM = ERRQ - TX2
!              IF (TEM < ZERO .AND. ERRQ > 0.1_kp) THEN
               IF (TEM < ZERO .AND. ERRQ > 0.5_kp) THEN
!              IF (TEM < ZERO .and.                                    &
!    &            (ntla < numtla .or. ERRQ > 0.5_kp)) THEN
                 SKPUP = .TRUE.                           ! No convergence !
                 ERRQ  = 10.0_kp                          ! No rain profile!
!!!!           ELSEIF (ABS(TEM) < ERRMI2 .OR. TX2 < ERRMIN) THEN
               ELSEIF (TX2 < ERRMIN) THEN
                 SKPUP = .TRUE.                           ! Converges      !
                 ERRQ = zero                              ! Rain profile exists!
               elseif (tem < zero .and. errq < 0.1_kp) then
                 skpup = .true.
!                if (ntla == numtla .or. tem > -0.003) then
                   errq  = zero
!                else
!                  errq = 10.0
!                endif
               ELSE
                 ERRQ = TX2                               ! Further iteration !

!              if (itr == itrmu .and. ERRQ > ERRMIN*10                  &
!    &            .and. ntla == 1) ERRQ = 10.0 
               ENDIF
            ENDIF
!
          ENDIF                                           ! SKPUP  ENDIF!
!
        ENDDO                                             ! End of the ITR Loop!!
!
        IF (ERRQ < 0.1_kp) THEN
          DDFT = .TRUE.
          RNB  = - RNB
!         do l=kd1,kb1-1
!           if (wvl(l)-wcbase < 1.0E-9) ddft = .false.
!         enddo
        ELSE
          DDFT = .FALSE.
        ENDIF

      enddo                                          ! End of ntla loop
!
!     Caution !! Below is an adjustment to rain flux to maintain
!                conservation of precip!
!
      IF (DDFT) THEN
        TX1 = zero
        DO L=KD,KB1
          TX1 = TX1 + RNF(L)
        ENDDO
        TX1 = TRAIN / (TX1+RNT+RNB)
        IF (ABS(TX1-one) < 0.2_kp) THEN
          RNT = MAX(RNT*TX1,ZERO)
          RNB = RNB * TX1
          DO L=KD,KB1
            RNF(L) = RNF(L) * TX1
          ENDDO
!      rain flux adjustment is over

        ELSE
          DDFT = .FALSE.
          ERRQ = 10.0_kp
        ENDIF
      ENDIF
!
      DOF = zero
      IF (.NOT. DDFT) then
        wvlu(kd:kp1) = zero
        RETURN     ! Rain profile did not converge!
                   ! No down draft for this case - rerurn
                   ! ------------------------------------
!
      else         ! rain profile converged - do downdraft calculation
                   ! ------------------------------------------------

        wvlu(kd:kp1) = wvl(kd:kp1) ! save updraft vertical velocity for output

!
!     Downdraft calculation begins
!     ----------------------------
!
        DO L=KD,K
          WCB(L) = zero
        ENDDO
!
        ERRQ  = 10.0_kp
! At this point stlt contains inverse of updraft vertical velocity 1/Wu.

        KK = MAX(KB1,KD1)
        DO L=KK,K
          STLT(L) = STLT(L-1)
        ENDDO
        TEM = stla / BB1       ! this is 2/(pi*radius*grav)
!
        DO L=KD,K
          IF (L <= KBL) THEN
            STLT(L) = ETA(L) * STLT(L) * TEM / ROR(L)
          ELSE
            STLT(L) = zero
          ENDIF
        ENDDO

        rsum1 = zero
        rsum2 = zero
!
        IDN(:) = idnmax
        DO L=KD,KP1
          ETD(L)  = zero
          WVL(L)  = zero
!         QRP(L)  = zero
        ENDDO
        DO L=KD,K
          EVP(L)   = zero
          BUY(L)   = zero
          QRP(L+1) = zero
        ENDDO
        HOD(KD)  = HOL(KD)
        QOD(KD)  = QOL(KD)
        TX1      = zero
!!!     TX1      = STLT(KD)*QRB(KD)*ONE              ! sigma at the top
!       TX1      = MIN(STLT(KD)*QRB(KD)*ONE, ONE)    ! sigma at the top
!       TX1      = MIN(STLT(KD)*QRB(KD)*0.5, ONE)    ! sigma at the top
        RNTP     = zero
        TX5      = TX1
        QA(1)    = zero
!
!       Here we assume RPART of detrained rain RNT goes to Pd
!
        IF (RNT > zero) THEN
          if (TX1 > zero) THEN
            QRP(KD) = (RPART*RNT / (ROR(KD)*TX1*GMS(KD)))               &
     &                                          ** (one/1.1364_kp)
          else
            tx1 = RPART*RNT / (ROR(KD)*GMS(KD)*QRP(KD)**1.1364_kp)
          endif
          RNTP    = (one - RPART) * RNT
          BUY(KD) = - ROR(KD) * TX1 * QRP(KD)
        ELSE
          QRP(KD) = zero
        ENDIF
!
!     L-loop for the downdraft iteration from KD1 to KP1 (bottom surface)
!
!       BUD(KD) = ROR(KD)
        idnm = 1
        DO L=KD1,KP1

          QA(1) = zero
          ddlgk = idn(idnm) == idnmax
          if (.not. ddlgk) cycle
          IF (L <= K) THEN
            ST1   = one - ALFIND(L)
            WA(1) = ALFIND(L)*HOL(L-1) + ST1*HOL(L)
            WA(2) = ALFIND(L)*QOL(L-1) + ST1*QOL(L)
            WA(3) = ALFIND(L)*TOL(L-1) + ST1*TOL(L)
            QA(2) = ALFIND(L)*HST(L-1) + ST1*HST(L)
            QA(3) = ALFIND(L)*QST(L-1) + ST1*QST(L)
          ELSE
            WA(1) = HOL(K)
            WA(2) = QOL(K)
            WA(3) = TOL(K)
            QA(2) = HST(K)
            QA(3) = QST(K)
          ENDIF
!
          FAC = two
          IF (L == KD1) FAC = one

          FACG    = FAC * half * GMF5     !  12/17/97
!
!         DDLGK   =  IDN(idnm) == 99

          BUD(KD) = ROR(L)

          TX1    = TX5
          WVL(L) = MAX(WVL(L-1),ONE_M1)

          QRP(L) = MAX(QRP(L-1),QRP(L))
!
!         VT(1)  = GMS(L-1) * QRP(L-1) ** 0.1364
          VT(1)  = GMS(L-1) * QRPF(QRP(L-1))
          RNT    = ROR(L-1) * (WVL(L-1)+VT(1))*QRP(L-1)
!
!         TEM    = MAX(ALM, 2.5E-4) * MAX(ETA(L), 1.0)
          TEM    = MAX(ALM,ONE_M6) * MAX(ETA(L), ONE)
!         TEM    = MAX(ALM, 1.0E-5) * MAX(ETA(L), 1.0)
          TRW(1) = PICON*TEM*(QRB(L-1)+QRT(L-1))
          TRW(2) = one / TRW(1)
!
          VRW(1) = half * (GAM(L-1) + GAM(L))
          VRW(2) = one / (VRW(1) + VRW(1))
!
          TX4    =  (QRT(L-1)+QRB(L-1))*(ONEBG*FAC*500.0_kp*EKNOB)
!
          DOFW   = one / (WA(3) * (one + NU*WA(2)))      !  1.0 / TVbar!
!
          ETD(L) = ETD(L-1)
          HOD(L) = HOD(L-1)
          QOD(L) = QOD(L-1)
!
          ERRQ   = 10.0_kp

!
          IF (L <= KBL) THEN
            TX3 = STLT(L-1) * QRT(L-1) * (half*FAC)
            TX8 = STLT(L)   * QRB(L-1) * (half*FAC)
            TX9 = TX8 + TX3
          ELSE
            TX3 = zero
            TX8 = zero
            TX9 = zero
          ENDIF
!
          TEM  = WVL(L-1) + VT(1)
          IF (TEM > zero) THEN
            TEM1 = one / (TEM*ROR(L-1))
            TX3 = VT(1) * TEM1 * ROR(L-1) * TX3
            TX6 = TX1 * TEM1
          ELSE
            TX6 = one
          ENDIF
!
          IF (L == KD1) THEN
            IF (RNT > zero) THEN
              TEM    = MAX(QRP(L-1),QRP(L))
              WVL(L) = TX1 * TEM * QRB(L-1)*(FACG*5.0_kp)
            ENDIF
            WVL(L) = MAX(ONE_M2, WVL(L))
            TRW(1) = TRW(1) * half
            TRW(2) = TRW(2) + TRW(2)
          ELSE
            IF (DDLGK) EVP(L-1) = EVP(L-2)
          ENDIF
!
!       No downdraft above level IDH
!

          IF (L < IDH) THEN

            ETD(L)   = zero
            HOD(L)   = WA(1)
            QOD(L)   = WA(2)
            EVP(L-1) = zero
            WVL(L)   = zero
            QRP(L)   = zero
            BUY(L)   = zero
            TX5      = TX9
            ERRQ     = zero
            RNTP     = RNTP + RNT * TX1
            RNT      = zero
            WCB(L-1) = zero

!         ENDIF
!         BUD(KD) = ROR(L)
!
!       Iteration loop for a given level L begins
!
          else
            DO ITR=1,ITRMD
!
!             cnvflg =  DDLGK .AND. (ERRQ > ERRMIN)
              cnvflg =  ERRQ > ERRMIN
              IF (cnvflg) THEN
!
!               VT(1) = GMS(L) * QRP(L) ** 0.1364
                VT(1) = GMS(L) * QRPF(QRP(L))
                TEM   =  WVL(L) + VT(1)
!
                IF (TEM > zero) THEN
                  ST1  = ROR(L) * TEM * QRP(L) + RNT
                  IF (ST1 /= zero) ST1 = two * EVP(L-1) / ST1
                  TEM1 = one / (TEM*ROR(L))
                  TEM2 = VT(1) * TEM1 * ROR(L) * TX8
                ELSE
                  TEM1 = zero
                  TEM2 = TX8
                  ST1  = zero
                ENDIF
!
                st2 = tx5
                TEM = ROR(L)*WVL(L) - ROR(L-1)*WVL(L-1)
                if (tem > zero) then
                  TX5 = (TX1 - ST1 + TEM2 + TX3)/(one+tem*tem1)
                else
                  TX5 = TX1 - tem*tx6 - ST1 + TEM2 + TX3
                endif
                TX5   = MAX(TX5,ZERO)
                tx5   = half * (tx5 + st2)
!
!               qqq = 1.0 + tem * tem1 * (1.0 - sialf)
!
!               if (qqq > 0.0) then
!                 TX5   = (TX1 - sialf*tem*tx6 - ST1 + TEM2 + TX3) / qqq
!               else
!                 TX5   = (TX1 - tem*tx6 - ST1 + TEM2 + TX3)
!               endif
!
                TEM1   = ETD(L)
                ETD(L) = ROR(L) * TX5 * MAX(WVL(L),ZERO)
!
                if (etd(l) > zero) etd(l) = half * (etd(l) + tem1)
!

                DEL_ETA = ETD(L) - ETD(L-1)

!                 TEM       = DEL_ETA * TRW(2)
!                 TEM2      = MAX(MIN(TEM, 1.0), -1.0)
!                 IF (ABS(TEM) > 1.0 .AND. ETD(L) > 0.0 ) THEN
!                   DEL_ETA = TEM2 * TRW(1)
!                   ETD(L)  = ETD(L-1) + DEL_ETA
!                 ENDIF
!                 IF (WVL(L) > 0.0) TX5 = ETD(L) / (ROR(L)*WVL(L))
!
                  ERRE  = ETD(L) - TEM1
!
                  tem  = max(abs(del_eta), trw(1))
                  tem2 = del_eta / tem
                  TEM1 = SQRT(MAX((tem+DEL_ETA)*(tem-DEL_ETA),ZERO))
!                 TEM1 = SQRT(MAX((TRW(1)+DEL_ETA)*(TRW(1)-DEL_ETA),0.0))

                  EDZ  = (half + ASIN(TEM2)*PIINV)*DEL_ETA + TEM1*PIINV

                DDZ   = EDZ - DEL_ETA
                WCB(L-1) = ETD(L) + DDZ
!
                TEM1  = HOD(L)
                IF (DEL_ETA > zero) THEN
                  QQQ    = one / (ETD(L) + DDZ)
                  HOD(L) = (ETD(L-1)*HOD(L-1) + DEL_ETA*HOL(L-1)        &
     &                                              + DDZ*WA(1)) * QQQ
                  QOD(L) = (ETD(L-1)*QOD(L-1) + DEL_ETA*QOL(L-1)        &
     &                                              + DDZ*WA(2)) * QQQ
                ELSEif((ETD(L-1) + EDZ) > zero) then
                  QQQ    = one / (ETD(L-1) + EDZ)
                  HOD(L) = (ETD(L-1)*HOD(L-1) + EDZ*WA(1)) * QQQ
                  QOD(L) = (ETD(L-1)*QOD(L-1) + EDZ*WA(2)) * QQQ
                ENDIF
                ERRH  = HOD(L) - TEM1
                ERRQ  = ABS(ERRH/HOD(L))  + ABS(ERRE/MAX(ETD(L),ONE_M5))

                DOF   = DDZ
                VT(2) = QQQ
!
                DDZ  = DOF
                TEM4 = QOD(L)
                TEM1 = VRW(1)
!
                QHS  = QA(3) + half * (GAF(L-1)+GAF(L)) * (HOD(L)-QA(2))
!
!                                           First iteration       !
!
                ST2  = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
                TEM2 = ROR(L) * QRP(L)
                CALL QRABF(TEM2,QRAF,QRBF)
                TEM6 = TX5 * (1.6_kp + 124.9_kp * QRAF) * QRBF * TX4
!
                CE   = TEM6 * ST2 / ((5.4e5_kp*ST2 + 2.55e6_kp)*(ETD(L)+DDZ))
!
                TEM2   = - ((one+TEM1)*(QHS+CE) + TEM1*QOD(L))
                TEM3   = (one + TEM1) * QHS * (QOD(L)+CE)
                TEM    = MAX(TEM2*TEM2 - four*TEM1*TEM3,ZERO)
                QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
!
!
!                                            second iteration   !
!
                ST2  = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
                CE   = TEM6 * ST2 / ((5.4e5_kp*ST2 + 2.55e6_kp)*(ETD(L)+DDZ))
!               CEE  = CE * (ETD(L)+DDZ)
!


                TEM2   = - ((one+TEM1)*(QHS+CE) + TEM1*tem4)
                TEM3   = (one + TEM1) * QHS * (tem4+CE)
                TEM    = MAX(TEM2*TEM2 - four*TEM1*TEM3,ZERO)
                QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
!                                              Evaporation in Layer L-1
!
                EVP(L-1) = (QOD(L)-TEM4) * (ETD(L)+DDZ)
!                                              Calculate Pd (L+1/2)
                QA(1)    = TX1*RNT + RNF(L-1) - EVP(L-1)
!
                if (qa(1) > zero) then
                  IF (ETD(L) > zero) THEN
                    TEM    = QA(1) / (ETD(L)+ROR(L)*TX5*VT(1))
                    QRP(L) = MAX(TEM,ZERO)
                  ELSEIF (TX5 > zero) THEN
                    QRP(L) = (MAX(ZERO,QA(1)/(ROR(L)*TX5*GMS(L))))      &
     &                                            ** (one/1.1364_kp)
                  ELSE
                    QRP(L) = zero
                  ENDIF
                else
                  qrp(l) = half * qrp(l)
                endif
!                                              Compute Buoyancy
                TEM1   = WA(3) + (HOD(L)-WA(1)-ALHL*(QOD(L)-WA(2)))     &
     &                         * onebcp

                TEM1   = TEM1 * (one + NU*QOD(L))
                ROR(L) = CMPOR * PRL(L) / TEM1
                TEM1   = TEM1 * DOFW
!!!             TEM1   = TEM1 * (1.0 + NU*QOD(L)) * DOFW

                BUY(L) = (TEM1 - one - QRP(L)) * ROR(L) * TX5
!                                              Compute W (L+1/2)

                TEM1   = WVL(L)
                WVL(L) = VT(2) * (ETD(L-1)*WVL(L-1) - FACG              &
     &                   * (BUY(L-1)*QRT(L-1)+BUY(L)*QRB(L-1)))

!
                if (wvl(l) < zero) then
!                 WVL(L) = max(wvl(l), 0.1*tem1)
!                 WVL(L) = 0.5*tem1
!                 WVL(L) = 0.1*tem1
!                 WVL(L) = 0.0
                  WVL(L) = 1.0e-10_kp
                else
                  WVL(L) = half*(WVL(L)+TEM1)
                endif

!
!               WVL(L) = max(0.5*(WVL(L)+TEM1), 0.0)

                ERRW   = WVL(L) - TEM1
!
                ERRQ   = ERRQ + ABS(ERRW/MAX(WVL(L),ONE_M5))


!               IF (ITR >= MIN(ITRMIN,ITRMD/2)) THEN
                IF (ITR >= MIN(ITRMND,ITRMD/2)) THEN

                  IF (ETD(L-1) == zero .AND. ERRQ > 0.2_kp) THEN

                    ROR(L) = BUD(KD)
                    ETD(L) = zero
                    WVL(L) = zero
                    ERRQ   = zero
                    HOD(L) = WA(1)
                    QOD(L) = WA(2)
!                   TX5      = TX1 + TX9
                    if (L <= KBL) then
                      TX5 = TX9
                    else
                      TX5 = (STLT(KB1) * QRT(KB1)                       &
     &                    +  STLT(KBL) * QRB(KB1)) * (0.5_kp*FAC)
                    endif

                    EVP(L-1) = zero
                    TEM      = MAX(TX1*RNT+RNF(L-1),ZERO)
                    QA(1)    = TEM - EVP(L-1)
!                   IF (QA(1) > 0.0) THEN

                    QRP(L)   = (QA(1) / (ROR(L)*TX5*GMS(L)))            &
     &                                              ** (one/1.1364_kp)
!                   endif
                    BUY(L)   = - ROR(L) * TX5 * QRP(L)
                    WCB(L-1) = zero
                  ENDIF
!
                  DEL_ETA = ETD(L) - ETD(L-1)
                  IF(DEL_ETA < zero .AND. ERRQ > 0.1_kp) THEN
                    ROR(L)   = BUD(KD)
                    ETD(L)   = zero
                    WVL(L)   = zero
!!!!!               TX5      = TX1 + TX9
                    CLDFRD(L-1) = TX5
!
                    DEL_ETA  = - ETD(L-1)
                    EDZ      = zero
                    DDZ      = -DEL_ETA
                    WCB(L-1) = DDZ
!
                    HOD(L)   = HOD(L-1)
                    QOD(L)   = QOD(L-1)
!
                    TEM4     = QOD(L)
                    TEM1     = VRW(1)
!
                    QHS      = QA(3) + half * (GAF(L-1)+GAF(L))         &
     &                                      * (HOD(L)-QA(2))

!
!                                           First iteration       !
!
                    ST2  = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
                    TEM2 = ROR(L) * QRP(L-1)
                    CALL QRABF(TEM2,QRAF,QRBF)
                    TEM6 = TX5 * (1.6_kp + 124.9_kp * QRAF) * QRBF * TX4
!
                    CE   = TEM6*ST2/((5.4e5_kp*ST2 + 2.55e6_kp)*(ETD(L)+DDZ))
!

                    TEM2   = - ((one+TEM1)*(QHS+CE) + TEM1*QOD(L))
                    TEM3   = (one + TEM1) * QHS * (QOD(L)+CE)
                    TEM    = MAX(TEM2*TEM2 -FOUR*TEM1*TEM3,ZERO)
                    QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))
!
!                                            second iteration   !
!
                    ST2  = PRL(L) * (QHS + TEM1 * (QHS-QOD(L)))
                    CE   = TEM6*ST2/((5.4e5_kp*ST2 + 2.55e6_kp)*(ETD(L)+DDZ))
!                   CEE  = CE * (ETD(L)+DDZ)
!


                    TEM2   = - ((one+TEM1)*(QHS+CE) + TEM1*tem4)
                    TEM3   = (one + TEM1) * QHS * (tem4+CE)
                    TEM    = MAX(TEM2*TEM2 -FOUR*TEM1*TEM3,ZERO)
                    QOD(L) = MAX(TEM4, (- TEM2 - SQRT(TEM)) * VRW(2))

!                                              Evaporation in Layer L-1
!
                    EVP(L-1) = (QOD(L)-TEM4) * (ETD(L)+DDZ)

!                                               Calculate Pd (L+1/2)
!                   RNN(L-1) = TX1*RNT + RNF(L-1) - EVP(L-1)

                    QA(1)    = TX1*RNT + RNF(L-1)
                    EVP(L-1) = min(EVP(L-1), QA(1))
                    QA(1)    = QA(1) - EVP(L-1)
                    qrp(l)   = zero

!
!                   IF (QA(1) > 0.0) THEN
!!                    RNS(L-1) = QA(1)
!!!                   tx5      = tx9
!                     QRP(L) = (QA(1) / (ROR(L)*TX5*GMS(L)))              &
!    &                                           ** (1.0/1.1364)
!                   endif
!                   ERRQ   = 0.0
!                                              Compute Buoyancy
!                   TEM1   = WA(3)+(HOD(L)-WA(1)-ALHL*(QOD(L)-WA(2)))     &
!    &                                                    * (1.0/CP)
!                   TEM1   = TEM1 * (1.0 + NU*QOD(L)) * DOFW
!                   BUY(L) = (TEM1 - 1.0 - QRP(L)) * ROR(L) * TX5
!
!                   IF (QA(1) > 0.0) RNS(L) = QA(1)

                    IF (L .LE. K) THEN
                       RNS(L) = QA(1)
                       QA(1)  = zero
                    ENDIF
                    tx5    = tx9
                    ERRQ   = zero
                    QRP(L) = zero
                    BUY(L) = zero
!
                  ENDIF
                ENDIF
              ENDIF
!
            ENDDO                ! End of the iteration loop  for a given L!
            IF (L <= K) THEN
              IF (ETD(L-1) == zero .AND. ERRQ > 0.1_kp .and. l <= kbl) THEN
!!!  &           .AND. ERRQ > ERRMIN*10.0 .and. l <= kbl) THEN
!    &           .AND. ERRQ > ERRMIN*10.0) THEN
                 ROR(L)   = BUD(KD)
                 HOD(L)   = WA(1)
                 QOD(L)   = WA(2)
                 TX5      =       TX9     ! Does not make too much difference!
!                TX5      = TX1 + TX9
                 EVP(L-1) = zero
!                EVP(L-1) = CEE * (1.0 - qod(l)/qa(3))
                 QA(1)    = TX1*RNT + RNF(L-1)
                 EVP(L-1) = min(EVP(L-1), QA(1))
                 QA(1)    = QA(1) - EVP(L-1)

!                QRP(L)   = 0.0
!              if (tx5 == 0.0 .or. gms(l) == 0.0) then
!                write(0,*)' Ctx5=',tx5,' gms=',gms(l),' ror=',ror(l)   &
!    &,          ' L=',L,' QA=',QA(1),' tx1=',tx1,' tx9=',tx9           &
!    &,          ' kbl=',kbl,' etd1=',etd(l-1),' DEL_ETA=',DEL_ETA
!              endif
!                IF (QA(1) > 0.0) THEN

                   QRP(L) = (QA(1) / (ROR(L)*TX5*GMS(L)))               &
     &                                           ** (one/1.1364_kp)
!                ENDIF
                 ETD(L)   = zero
                 WVL(L)   = zero
                 ST1      = one - ALFIND(L)

                 ERRQ     = zero
                 BUY(L)   = - ROR(L) * TX5 * QRP(L)
                 WCB(L-1) = zero
              ENDIF
            ENDIF
!
            LL = MIN(IDN(idnm), KP1)
            IF (ERRQ < one .AND. L <= LL) THEN
              IF (ETD(L-1) > zero .AND. ETD(L) == zero) THEN
               IDN(idnm) = L
               wvl(l)    = zero
               if (L < KBL .or. tx5 > zero) idnm  = idnm + 1
               errq = zero
              ENDIF
              if (etd(l) == zero .and. l > kbl) then
                idn(idnm) = l
                if (tx5 > zero) idnm  = idnm + 1
              endif
            ENDIF
! 
!     If downdraft properties are not obtainable, (i.e.solution does
!      not converge) , no downdraft is assumed
!
!           IF (ERRQ > ERRMIN*100.0 .AND. IDN(idnm) == 99)                &
            IF (ERRQ > 0.1_kp .AND. IDN(idnm) == idnmax) DDFT = .FALSE.
!
            DOF = zero
            IF (.NOT. DDFT) RETURN
!
!         if (ddlgk .or. l .le. idn(idnm)) then
!           rsum2 = rsum2 + evp(l-1)
!           write(0,*)' rsum1=',rsum1,' rsum2=',rsum2,' L=',L,' qa=',qa(1)&
!    &,   ' evp=',evp(l-1)
!         else
!           rsum1 = rsum1 + rnf(l-1)
!           write(0,*)' rsum1=',rsum1,' rsum2=',rsum2,' L=',L,' rnf=',    &
!     &     rnf(l-1)
!         endif

          endif         ! if (l < idh)
        ENDDO                      ! End of the L Loop of downdraft !

        TX1 = zero

        DOF = QA(1)
!
!       write(0,*)' dof=',dof,' rntp=',rntp,' rnb=',rnb
!       write(0,*)' total=',(rsum1+dof+rntp+rnb)
!
        dof     = max(dof, zero)
        RNN(KD) = RNTP
        TX1     = EVP(KD)
        TX2     = RNTP + RNB + DOF

        II = IDH
        IF (II >= KD1+1) THEN
           RNN(KD)   = RNN(KD) + RNF(KD)
           TX2       = TX2 + RNF(KD)
           RNN(II-1) = zero
           TX1       = EVP(II-1)
        ENDIF
        DO L=KD,K
          II = IDH

          IF (L > KD1 .AND. L < II) THEN
            RNN(L-1) = RNF(L-1)
            TX2      = TX2 + RNN(L-1)
          ELSEIF (L >= II .AND. L < IDN(idnm)) THEN
            rnn(l)   = rns(l)
            tx2      = tx2 + rnn(l)
            TX1      = TX1 + EVP(L)
          ELSEIF (L >= IDN(idnm)) THEN
            ETD(L+1) = zero
            HOD(L+1) = zero
            QOD(L+1) = zero
            EVP(L)   = zero
            RNN(L)   = RNF(L) + RNS(L)
            TX2      = TX2    + RNN(L)
          ENDIF
        ENDDO
!
!      For Downdraft case the rain is that falls thru the bottom

        L = KBL

        RNN(L)    = RNN(L) + RNB
        CLDFRD(L) = TX5

!
!     Caution !! Below is an adjustment to rain flux to maintain
!                conservation of precip!

!
        IF (TX1 > zero) THEN
          TX1 = (TRAIN - TX2) / TX1
        ELSE
          TX1 = zero
        ENDIF

        DO L=KD,K
          EVP(L) = EVP(L) * TX1
        ENDDO

      ENDIF                       ! if (.not. DDFT) loop   endif
!
!***********************************************************************
!***********************************************************************

      RETURN
      end subroutine ddrft

!>\ingroup rascnv_schm
      SUBROUTINE QSATCN(TT,P,Q,DQDT)
!
      USE FUNCPHYS , ONLY : fpvs

      implicit none
!
      real(kind=kind_phys) TT, P, Q, DQDT
!
!     real(kind=kind_phys), parameter :: ZERO=0.0, ONE=1.0              &
!    &,                                  rvi=one/rv,     facw=CVAP-CLIQ &
!    &,                                  faci=CVAP-CSOL, hsub=alhl+alhf &
!    &,                                  tmix=TTP-20.0                  &
!    &,                                  DEN=one/(TTP-TMIX)
!
      real(kind=kind_phys) es, d, hlorv, W
!
!     es = 10.0 * fpvs(tt)                ! fpvs is in centibars!
      es = min(p, 0.01_kp * fpvs(tt))     ! fpvs is in Pascals!
!     D  = one / max(p+epsm1*es,ONE_M10)
      D  = one / (p+epsm1*es)
!
      q  = MIN(eps*es*D, ONE)
!
      W     = max(ZERO, min(ONE, (TT - TMIX)*DEN))
      hlorv = ( W      * (alhl + FACW * (tt-ttp))                       &
     &       + (one-W) * (hsub + FACI * (tt-ttp)) ) * RVI
      dqdt  = p * q * hlorv *  D / (tt*tt)
!
      return
      end subroutine qsatcn

!>\ingroup rascnv_schm
      SUBROUTINE ANGRAD(PRES, ALM,  AL2, TLA)
      implicit none

      real(kind=kind_phys) PRES, ALM,  AL2,  TLA,  TEM
!
      integer i
!
      IF (TLA < 0.0_kp) THEN
          IF (PRES <= PLAC(1)) THEN
            TLA = TLAC(1)
          ELSEIF (PRES <= PLAC(2)) THEN
            TLA = TLAC(2) + (PRES-PLAC(2))*tlbpl(1)
          ELSEIF (PRES <= PLAC(3)) THEN
            TLA = TLAC(3) + (PRES-PLAC(3))*tlbpl(2)
          ELSEIF (PRES <= PLAC(4)) THEN
            TLA = TLAC(4) + (PRES-PLAC(4))*tlbpl(3)
          ELSEIF (PRES <= PLAC(5)) THEN
            TLA = TLAC(5) + (PRES-PLAC(5))*tlbpl(4)
          ELSEIF (PRES <= PLAC(6)) THEN
            TLA = TLAC(6) + (PRES-PLAC(6))*tlbpl(5)
          ELSEIF (PRES <= PLAC(7)) THEN
            TLA = TLAC(7) + (PRES-PLAC(7))*tlbpl(6)
          ELSEIF (PRES <= PLAC(8)) THEN
            TLA = TLAC(8) + (PRES-PLAC(8))*tlbpl(7)
          ELSE
            TLA = TLAC(8)
          ENDIF
      ENDIF
        IF (PRES >= REFP(1)) THEN
          TEM = REFR(1)
        ELSEIF (PRES >= REFP(2)) THEN
          TEM = REFR(1) + (PRES-REFP(1)) * drdp(1)
        ELSEIF (PRES >= REFP(3)) THEN
          TEM = REFR(2) + (PRES-REFP(2)) * drdp(2)
        ELSEIF (PRES >= REFP(4)) THEN
          TEM = REFR(3) + (PRES-REFP(3)) * drdp(3)
        ELSEIF (PRES >= REFP(5)) THEN
          TEM = REFR(4) + (PRES-REFP(4)) * drdp(4)
        ELSEIF (PRES >= REFP(6)) THEN
          TEM = REFR(5) + (PRES-REFP(5)) * drdp(5)
        ELSE
          TEM = REFR(6)
        ENDIF
!
        tem = 2.0e-4_kp / tem
        al2 = min(4.0_kp*tem, max(alm, tem))
!
      RETURN
      end subroutine angrad

!>\ingroup rascnv_schm
      SUBROUTINE SETQRP
      implicit none

      real(kind=kind_phys) tem2,tem1,x,xinc,xmax,xmin
      integer jx
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     XMIN   = 1.0E-6
      XMIN   = 0.0_kp
      XMAX   = 5.0_kp
      XINC   = (XMAX-XMIN)/(NQRP-1)
      C2XQRP = one / XINC
      C1XQRP = one - XMIN*C2XQRP
      TEM1   = 0.001_kp ** 0.2046_kp
      TEM2   = 0.001_kp ** 0.525_kp
      DO JX=1,NQRP
        X         = XMIN + (JX-1)*XINC
        TBQRP(JX) =        X ** 0.1364_kp
        TBQRA(JX) = TEM1 * X ** 0.2046_kp
        TBQRB(JX) = TEM2 * X ** 0.525_kp
      ENDDO    
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      end subroutine setqrp

!>\ingroup rascnv_schm
      SUBROUTINE QRABF(QRP,QRAF,QRBF)
      implicit none
!
      real(kind=kind_phys) QRP, QRAF, QRBF, XJ, REAL_NQRP
      INTEGER JX
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      REAL_NQRP = REAL(NQRP)
      XJ   = MIN(MAX(C1XQRP+C2XQRP*QRP,ONE),REAL_NQRP)
      JX   = MIN(XJ,NQRP-ONE)
      XJ   = XJ - JX
      QRAF = TBQRA(JX)  + XJ * (TBQRA(JX+1)-TBQRA(JX))
      QRBF = TBQRB(JX)  + XJ * (TBQRB(JX+1)-TBQRB(JX))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      end subroutine qrabf

!>\ingroup rascnv_schm
      SUBROUTINE SETVTP
      implicit none

      real(kind=kind_phys), parameter :: vtpexp=-0.3636_kp
      real(kind=kind_phys) xinc,x,xmax,xmin
      integer jx
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      XMIN   = 0.05_kp
      XMAX   = 1.5_kp
      XINC   = (XMAX-XMIN)/(NVTP-1)
      C2XVTP = one / XINC
      C1XVTP = one - XMIN*C2XVTP
      DO JX=1,NVTP
        X         = XMIN + (JX-1)*XINC
        TBVTP(JX) =        X ** VTPEXP
      ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      end subroutine setvtp
!
!>\ingroup rascnv_schm
      real(kind=kind_phys) FUNCTION QRPF(QRP)
!
      implicit none

      real(kind=kind_phys) QRP, XJ, REAL_NQRP
      INTEGER JX
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      REAL_NQRP = REAL(NQRP)
      XJ   = MIN(MAX(C1XQRP+C2XQRP*QRP,ONE),REAL_NQRP)
!     XJ   = MIN(MAX(C1XQRP+C2XQRP*QRP,ONE),FLOAT(NQRP))
      JX   = MIN(XJ,NQRP-ONE)
      QRPF = TBQRP(JX)  + (XJ-JX) * (TBQRP(JX+1)-TBQRP(JX))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      end function qrpf

!>\ingroup rascnv_schm
      real(kind=kind_phys) FUNCTION VTPF(ROR)
!
      implicit none
      real(kind=kind_phys) ROR, XJ, REAL_NVTP
      INTEGER JX
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      REAL_NVTP = REAL(NVTP)
      XJ   = MIN(MAX(C1XVTP+C2XVTP*ROR,ONE),REAL_NVTP)
      JX   = MIN(XJ,NVTP-ONE)
      VTPF = TBVTP(JX)  + (XJ-JX) * (TBVTP(JX+1)-TBVTP(JX))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      end function vtpf

!>\ingroup rascnv_schm
      real(kind=kind_phys) FUNCTION CLF(PRATE)
!
      implicit none
      real(kind=kind_phys) PRATE
!
      real (kind=kind_phys), parameter :: ccf1=0.30_kp, ccf2=0.09_kp    &
     &,                                   ccf3=0.04_kp, ccf4=0.01_kp    &
     &,                                   pr1=1.0_kp,   pr2=5.0_kp      &
     &,                                   pr3=20.0_kp
!
      if (prate < pr1) then
        clf = ccf1
      elseif (prate < pr2) then
        clf = ccf2
      elseif (prate < pr3) then
        clf = ccf3
      else
        clf = ccf4
      endif
!
      RETURN
      end function clf
      end module rascnv
